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Background: 

Considerable  research  has  shown  that  the  major  control  on  the  transport  and  fate  of  a 
pollutant  as  it  moves  through  an  aquifer  is  the  spatial  distribution  of  hydraulic 
conductivity.  Although  chemical  and  microbial  processes  play  important  roles,  their 
influence  cannot  be  understood  without  a  detailed  knowledge  of  the  subsurface  variations 
in  hydraulic  conductivity  at  a  site.  Many  theories  have  been  developed  to  quantify,  in  a 
generic  sense,  the  influence  of  these  variations  using  stochastic  processes  or  fractal 
representations.  It  is  increasingly  apparent,  however,  that  site-specific  features  of  the 
hydraulic  conductivity  distribution  (such  as  high  conductivity  zones)  need  to  be 
quantified  to  reliably  predict  contaminant  movement.  Conventional  hydraulic  field 
techniques  only  provide  information  of  a  highly  averaged  nature  or  information  restricted 
to  the  immediate  vicinity  of  the  test  well.  Therefore,  development  of  new  innovative 
methods  to  delineate  the  detailed  hydraulic  conductivity  distribution  at  a  given  site  should 
be  a  high  priority.  The  research  proposed  here  is  directed  at  addressing  this  problem  by 
developing  techniques  to  map  3-D  hydraulic  conductivity  distributions. 


Objective: 

Since  spatial  changes  in  hydraulic  conductivity  are  a  major  factor  governing  the  transport 
and  fate  of  a  pollutant  as  it  moves  through  an  aquifer,  we  focus  on  the  development  of 
new  innovative  methods  to  delineate  these  spatial  changes.  The  objective  of  the  research 
proposed  here  is  to  build  on  our  previous  work  to  develop  and  improve  field  techniques 
for  better  definition  of  the  three-dimensional  spatial  distribution  of  hydraulic  conductivity 
by  using  hydraulic  tomography  coupled  with  high-resolution  slug  testing. 


Technology  Approach: 

We  have  worked  for  many  years  to  quantify  hydraulic  conductivity  fields  in 
heterogeneous  aquifers.  One  promising  method  we  have  worked  on  extensively  is  high- 
resolution  slug  testing.  This  method  allows  the  delineation  of  the  vertical  distribution  of 
hydraulic  conductivity  near  an  observation  well.  We  propose  to  combine  this  method 
with  another  innovative  method  for  investigating  the  hydraulic  conductivity  distribution 
between  wells,  called  hydraulic  tomography.  We  will  use  an  oscillating  signal  and 
measure  its  phase  and  amplitude  through  space  in  order  to  estimate  the  hydraulic 
conductivity  distribution  of  the  material  through  which  it  has  traveled.  Our  preliminary 
work  shows  that  the  phase  and  amplitude  of  the  received  signal  can  be  measured  over 
reasonable  distances.  The  high-resolution  slug  testing  results  will  be  used  as  an  initial 
condition  and  will  provide  conditioning  for  the  tomographic  inverse  procedure,  to  help 
with  any  non-uniqueness  problems.  Slug  test  data  are  most  accurate  near  the  tested  well 
and  should  probably  not  be  extrapolated  blindly  between  wells.  Together,  slug  testing 
and  hydraulic  tomography  should  be  more  powerful  than  either  one  used  alone  and 
should  give  the  best  opportunity  to  characterize  the  hydraulic  conductivity  in-situ  by  a 
direct  measure  of  water  flow,  as  an  alternative  to  indirect  methods  using  geophysical 
techniques. 


4 


Introduction 


A  typical  method  used  to  determine  fluid  behavior  in  a  geologic  matrix  near  a 
well  is  a  pumping  test.  Here  a  pump  is  installed  into  a  well  and  groundwater  is  removed 
or  injected  while  water  levels  in  surrounding  observation  wells  are  monitored.  Then  the 
aquifer  parameters  can  be  estimated  by  monitoring  changes  in  water  levels  at  observation 
wells  at  some  distance.  These  tests  are  typically  large  in  scale  (Schad  and  Teutsch, 

1994).  Another  test  is  an  interference  test,  which  is  a  special  pumping  test  where  the 
pump  discharge  has  a  variable  rate.  Interference  tests  are  conducted  by  variable 
production  or  injection  of  fluid  (hydraulic  head  changes)  at  one  well,  and  observing  the 
changing  pressure  or  hydraulic  head  with  time  and  distance  at  other  locations.  These 
tests  are  valued  to  estimate  flow  characteristics  in  situ,  but  are  measures  of  the  aquifer 
material  over  large  volumes  also. 

On  the  other  hand,  physical  cores  of  aquifer  material  can  be  obtained  by  various 
drilling  methods.  These  samples  can  then  be  tested  in  a  laboratory  (i.e.,  falling  or 
constant  head  permeability  tests)  to  estimate  the  hydraulic  properties.  One  advantage  to 
this  method  is  that  the  sample  can  be  visually  inspected.  Some  disadvantages  to  this 
method  are  that  the  material  is  disturbed  from  its  natural  environment  and  the  sample  is  a 
small  representation  of  the  total  aquifer. 

Another  common  technique  for  determining  aquifer  parameters  is  slug  tests.  A 
slug  test  initiates  a  head  change  in  a  well,  then  monitors  the  response  of  the  aquifer 
material  to  estimate  the  hydraulic  conductivity  (K).  Slug  testing  is  usually  only 
conducted  in  a  single  well.  It  is  generally  accepted  that  the  radius  of  influence  of  a  slug 
test  is  small  and  only  provides  a  limited  view  of  subsurface  hydrogeologic  properties  near 
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the  well.  Traditionally,  slug  tests  have  been  initiated  with  the  addition  into  a  well  of  a 
known  volume  of  water  or  a  physical  slug.  More  recently,  pneumatic  methods  have 
become  popular  (Zemansky  and  McElwee,  2005;  Sellwood,  2001;  McCall  et  al.,  2000) 
for  multilevel  slug  testing.  Slug  tests  in  low  K  formations  can  take  much  longer  than  in 
material  with  high  permeability.  To  overcome  this,  the  fluid  column  in  a  well  can  be 
pressurized  and  the  pressure  change  with  time  can  be  used  as  an  alternative  (Bredehoeft 
and  Papadopulos,  1980). 


Figure  1 .  High  resolution  slug  testing  equipment  deployed 
in  a  fully  penetrating  well. 

Typical  slug  tests  are  conducted  by  exciting  the  entire  length  of  the  well  screen. 
Whole  well  slug  testing  can  provide  information  near  the  tested  well  but  it  is  an  average 
response  over  the  total  length  of  that  well’s  screen.  However,  aquifers  are  naturally 
heterogeneous  and  whole  well  slug  testing  is  unable  to  distinguish  areas  of  high  or  low  K. 
High  resolution  slug  testing  [(HRST),  over  short  screen  intervals  (Figure  1)],  provides  a 
more  detailed  vertical  profde  of  K  near  the  tested  well.  In  this  research  the  HRST 
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interval  is  approximately  0.5  m;  but,  stressed  intervals  as  small  as  5  cm  have  been  used 
(Healey  et  al.,  2004).  Currently  there  is  no  accepted  method  to  bridge  the  gap  between 
the  larger  lateral  well-to-well  averages  from  pumping  or  interference  tests  and  detailed 
vertical  estimates  of  K  from  HRST.  Proposed  here  is  a  method  to  obtain  estimates  of 
aquifer  parameters  at  larger  radii  of  influence,  while  simultaneously  maintaining  a  higher 
resolution. 

Pulse  testing  is  one  method  of  determining  fluid  flow  parameters  that  is  often 
employed  by  the  petroleum  industry.  Johnson  et  al.  (1966)  published  results  of 
experiments  conducted  in  a  sandstone  reservoir  near  Chandler,  OK.  They  found  that  the 
new  pulse  method  was  as  effective  as  typical  interference  tests.  The  transient  pressure 
signal  is  propagated  by  in  situ  fluid  and  is  therefore  a  direct  measure  of  reservoir 
diffusivity.  Other  advantages  of  the  pulse  method  are  the  ability  to  distinguish  the  test 
from  background  noise  because  of  its  controlled  frequency  of  oscillation  and  the 
reduction  of  down  time  relative  to  production.  Since  1966,  pulse  testing  has  been  used  to 
delineate  fractures  (Barker,  1988;  Brauchler,  et  al.,  2001)  and  to  predict  water  flood 
performance  (Pierce,  1977). 

Other  pulse  test  examples  include  tidal,  seismic,  and  oil  field  methods.  The 
changes  in  groundwater  levels  as  a  result  of  tidal  fluctuations  have  been  well  studied 
(Ferris,  1951;  Hantush,  1960;  Jiao  and  Tang,  1999).  The  sinusoidal  tidal  fluctuations  that 
propagate  inland  through  an  aquifer  are  related  to  aquifer  storativity  and  transmissivity. 
Solutions  to  water  level  fluctuations  induced  by  seismic  waves  were  presented  by  Cooper 
et  al.  (1965).  The  pressure  head  fluctuations  controlling  water  levels  are  a  result  of  the 
vertical  motion  of  the  aquifer  but  are  dominated  by  dilation  of  aquifer  porosity.  An 
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interference  test  of  alternating  oil  production  and  shut-in  time  was  conducted  to 
determine  the  interconnectivity  of  wells  in  a  production  field  (Johnson  et  ah,  1966).  Here 
the  source  well  is  assumed  to  be  a  line  source  in  an  infinite  homogeneous  reservoir.  The 
time  lag  and  the  received  amplitude  were  used  to  estimate  the  average  well-to-well 
transmissivity  and  storage  properties  of  the  reservoir.  These  oil  field  methods  were 
theoretically  adapted  to  hydrogeo  logic  characterization  by  Black  and  Kipp  (1981). 
Analytical  solutions  of  a  fracture  responding  to  a  single  pulse  interference  test,  a  slug  of 
water,  was  modeled  and  tested  by  Novakowski  (1989).  Straddle  packers  isolated  the 
fracture  and  were  used  to  apply  the  slug  of  water  by  being  deflated.  The  duration  of  these 
tests  was  on  average  30  min.  The  sequential  pumping  or  removal  of  water  was  used  to 
collect  head  responses  between  wells  (Yeh  and  Liu,  2000).  In  these  experiments  multiple 
ray  paths  were  analyzed  as  a  hydraulic  tomography  experiment.  Such  experiments  show 
promise  in  their  ability  to  distinguish  lateral  and  vertical  2-D  variations  in  heterogeneity 
by  changes  in  the  signal  over  the  travel  path. 

The  research  presented  in  this  report  uses  continuous,  controlled,  sinusoidal 
pressure  signals  [the  continuous  pulse  test  (CPT)]  as  a  means  to  estimate  vertical  profiles 
of  well-to-well  averaged  hydraulic  diffusivity.  In  this  research,  the  primary  method  of 
stimulation  of  the  alluvial  aquifer  was  achieved  by  pneumatic  methods  or  by  mechanical 
pumping  methods.  In  the  pneumatic  method  the  column  of  air  within  a  well  was 
pressurized  via  an  air  compressor.  A  signal  generator  or  computer  controlled  switch  was 
used  to  open  and  close  valves  at  the  well-head  allowing  air  to  enter  or  exit  the  well.  The 
signal  generator  or  computer  controlled  signal  produced  an  adjustable  frequency 
excitation  voltage,  which  controlled  the  periodicity  of  the  continuous  pulse-testing  signal. 
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Theoretically,  a  square  wave  pressure  test  is  the  simplest  to  conduct  because  of  the 
instantaneous  pressure  changes  (Lee,  1982).  Due  to  the  input  air  pressure,  the  water 
column  in  a  well  will  be  depressed  creating  flow  through  the  well  screen.  This  pulse  of 
hydraulic  pressure  is  transferred  to  the  aquifer  system  based  on  the  diffusivity  of  the 
material.  As  the  air  column  within  the  well  is  allowed  to  return  to  atmospheric  pressure, 
water  rushes  back  into  the  well  from  the  aquifer.  The  mechanical  pumping  method  used 
a  surface  reservoir  of  water  and  a  pump  to  inject  water  into  the  aquifer.  The  pressurized 
water  was  allowed  to  flow  into  the  aquifer  in  a  periodic  sinusoidal  fashion  with  the  help 
of  a  computer  controlled  valve.  These  fluctuations  are  periodic  and  similar  to  tidal 
fluctuations  acting  upon  a  costal  aquifer  system.  The  governing  equations  for  an  aquifer 
responding  to  tidal  fluctuations  were  adapted  to  Cartesian,  cylindrical,  and  spherical 
coordinate  systems  describing  groundwater  flow  with  sinusoidal  boundary  conditions,  in 
order  to  describe  the  data  used  in  this  report. 

The  period,  the  phase,  and  the  amplitude  of  the  produced  wave  can  then  be 
measured  simultaneously  at  the  source  well  and  at  observation  wells.  Through 
dispersion,  the  aquifer  material  will  decrease  the  fidelity  of  the  input  signal,  retard  the 
propagation,  and  attenuate  the  propagating  wave  front,  resulting  in  a  phase  lag  or  shift, 
and  a  decrease  in  the  amplitude.  The  amplitude  ratio  [received  amplitude  Ar  divided  by 
the  initial  amplitude  Ao]  and  the  phase  difference  [reference  phase  4>o  minus  the  received 
phase  (j)r]  can  then  be  used  to  calculate  the  hydraulic  diffusivity  (Lee,  1982). 

Zero  Offset  Profile  (ZOP,  source  and  receiver  at  same  elevation)  data  and 
Multiple  Offset  Gather  (MOG,  source  location  fixed;  receiver  elevation  varied)  data  were 
collected  at  the  University  of  Kansas’  Geohydro  logic  Experimental  and  Monitoring  Site 
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(GEMS),  a  well-studied  shallow  semi-confined  alluvial  aquifer  system  in  the  Kansas 
River  floodplain.  Line  sources  equal  to  the  total  screen  length  and  point  sources  isolated 
by  custom  bladder  packers  were  used  in  these  experiments.  Field  data  indicate  that 
sinusoidal  signals  can  propagate  reasonable  distances,  and  may  provide  estimates  of  the 
well-to-well  diffusivity.  Vertical  profiles  of  hydraulic  conductivity  (K),  measured  with 
high-resolution  slug  testing  (HRST),  were  collected  for  correlation  with  the  CPT  data. 

The  GEMS  area  is  located  in  Douglas  County,  northeast  Kansas,  along  the 
northern  margin  of  the  Kansas  River  floodplain  (Figure.  2).  GEMS  is  in  a  Pennsylvanian 
bedrock  valley  filled  with  Wisconsinan-age  glaciofluvial  terrace  sediments  (O’Conner, 
1960;  Schulmeister,  2000).  The  upper  1 1  m  of  sediments  are  mostly  silts  and  clays  and 
the  lower  12  m  of  sediments  at  GEMS  is  a  fining  upward  sequence  of  pebbles,  coarse 
sand,  and  fine  sand,  underlain  by  the  Tonganoxie  Sandstone  member  (Jiang,  1991). 
Within  the  sequences  of  sandy  material  are  lenses  of  low  permeability  fine-grained 
sediments.  These  clay  lenses  occur  at  various  elevations  and  can  be  up  to  1  m  thick 
(Schulmeister,  2000;  Healey  et  al.,  2004).  As  an  aquifer,  the  Kansas  River  alluvium  is  a 
prolific  deposit  of  unconsolidated  sands  and  gravels.  This  high  yielding  semi-confined 
aquifer  meets  the  needs  of  agricultural,  industrial,  and  community  interests. 

Many  studies  have  been  conducted  at  GEMS  and  many  well  nests  have  been 
completed  to  various  depths  with  various  screen  lengths.  Porosity,  grain  size,  and  K 
were  estimated  by  laboratory  experiments  performed  on  physical  samples  of  the  aquifer 
material  (Jiang,  1991).  A  single-well  injection  tracer  test  was  used  to  estimate  a  K 
distribution  by  monitoring  the  transport  of  an  electrolytic  solution  (Huettl,  1992).  The  K 
distribution  in  an  area  of  GEMS  was  also  estimated  by  conducting  an  induced-gradient 
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tracer  test  through  a  multilevel  groundwater  sampling  well  field  (Bohling,  1999).  Direct 
push  bulk  electrical  conductivity  (EC)  profiling  (Figure  3)  and  direct  push  pneumatic  slug 
tests  were  also  done  adjacent  to  the  tracer  experiment  well  field  (Sellwood,  2001). 


Figure  2.  GEMS  location  map  and  aerial  photographs. 
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Figure  3.  Direct  push  drilling  unit,  Electrical  Conductance  probe,  and  example  profile. 
Most  recently,  HRST  K  estimates  were  collected  in  numerous  wells  that  were  fully 
screened  through  the  aquifer  material  (Ross,  2004;  Ross  and  McElwee,  2007).  These 
independent  studies  and  the  research  presented  here  produced  estimates  of  K  that  can  be 
collected  into  a  database.  After  compiling  these  data,  vertical  and  lateral  variations  of  the 
K  distribution  are  evident.  Typically  at  GEMS,  K  increases  with  depth  in  the  sands  and 
gravels,  and  low  K  material  can  be  associated  with  high  EC  measurements,  usually 
associated  with  the  overlying  silt  and  clay  sediments.  In  most  areas  at  GEMS,  “layers”  or 
zones  of  high  K  material  are  apparent  in  the  sand  and  gravel  aquifer. 

Theory 

Fluid  flow  in  saturated  aquifers  behaves  much  like  heat  flow  and  can  be  described 
by  similar  equations.  Excess  pore  pressures,  matrix  permeability,  compressibility,  and 
storativity  all  influence  the  fluctuations  of  groundwater  levels  in  response  to  applied 
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stresses.  The  excess  fluid  pressure  Pe,  above  hydrostatic  pressure  Ps,  is  related  to  the  total 
stress  on  the  aquifer  g,  and  changes  the  stress  Ag  by 
(1)  G  +  AG  =  Ge  +  (Ps  +  Pe) 


The  above  equation  allocates  the  additional  stress  to  either  the  aquifer  matrix 
itself  (Ge )  or  to  excess  hydraulic  pressure,  Pe.  By  changing  the  hydraulic  pressure  or 
hydraulic  head,  the  water  levels  in  an  aquifer  also  change  accordingly.  The  total 
hydraulic  head  (h)  hydraulic  potential  measured  in  a  well  is  a  combination  of  the 
elevation  head  z,  and  the  hydraulic  pressure  head,  P 
(2)  h  =  z  +  P/pg 

such  that 


(3)  P  =  Ps  +  Pe 

Since  the  elevation  is  static,  the  only  dynamic  portion  of  h  is  due  to  pressure 
changes  as  shown  in  the  following  equation 


(4) 


dh  1  dP 
dt  pg  dt 


where  p  is  the  fluid  density  and  g  is  the  acceleration  of  gravity.  Substituting  equation  (3) 
into  equation  (2)  the  total  head  measured  in  a  well  can  also  be  expressed  as 
(5)  h  =  Z  +  (Ps/Pwg  +  Pe/pwg) 

Darcy’s  law  states  that  the  discharge  Q  of  a  fluid  through  a  porous  media  depends  on  the 

dh 

hydraulic  gradient  (the  change  in  head  with  distance)  — ,  and  the  cross  sectional  area  A. 

& L 


Darcy’s  Law  is 

(6) 
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The  above  equations  can  be  generalized  to  three  dimensions.  The  goal  of  this  research  is 
to  utilize  the  response  of  hydrogeologic  material  to  cyclic  pressure  signals  to  estimate  the 
D  or  K  distribution  in  an  aquifer. 

Groundwater  fluctuations  near  coastal  regions  have  been  studied  and  elementary 
equations  have  been  developed  to  associate  regional  groundwater  levels  with  tidal 
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fluctuations  (Hantush,  1960).  The  basic  mathematical  description  of  a  one-dimensional 
transient  pressure  head  signal  with  sinusoidal  boundary  conditions  [sin(27ift)]  is 


(12) 


h(r,t)  =  hoed  sin(<D0 -Or) . 


The  head  at  some  distance  and  time  h(r,t)  is  the  initial  amplitude  h0,  some  decay  term  ed, 
multiplied  by  the  sine  of  the  source  reference  phase  (O0=27ift)  minus  the  phase  shift,  <t>,. 
The  amplitude  decay  and  the  phase  shift  depend  on  the  ability  of  the  aquifer  to  transmit 
the  sinusoidal  signal.  Namely,  it  is  the  hydraulic  diffusivity  (D  or  K/Ss)  of  the  aquifer 
that  influences  the  hydraulic  head  measured  at  some  distance  and  time  from  the  source  of 
a  pressure  head  fluctuation.  Three  equations  for  the  head  response  to  the  propagation  of  a 
sinusoidal  boundary  condition  (causing  excess  fluid  pressure)  within  a  homogeneous 
isotropic  formation  have  been  adapted  from  equation  (12).  Equation  (12)  has  been 
extended  to  various  coordinate  systems,  presented  below. 

Linear  Cartesian  System 

(13)  h(x,t)  =  hae 

Cylindrical  Radial  System 


sin 


2  nft- 


K 


X 


(14) 


h(r,t)  =  h0 


( 


sin 


2nft  — 


V 


A 


J 


Spherical  Radial  System 


(15) 


h{r  ,t)  =  h0 


-sm 


27rft~i^tr 
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Where  t  is  time,  x  or  r  is  the  distance  from  the  source,  f  is  the  frequency,  h0  is  the  initial 
amplitude  of  the  pressure  head  fluctuation  at  the  source,  Ss  is  the  specific  storage,  and  K 
is  the  hydraulic  conductivity.  Specific  storage  is  the  volume  of  fluid  added  or  released 
per  unit  volume  of  aquifer  per  unit  thickness,  from  compression  or  relaxation  of  the 
aquifer  skeleton  and  pores  due  to  changes  in  stress.  Equation  (13)  is  exact  (Hantush, 
1960);  however,  equations  (14)  and  (15)  are  good  approximations  away  from  the  origin. 
This  issue  will  be  confirmed  by  numerical  modeling  later.  The  coordinate  equations  (13, 
14,  and  15)  can  be  thought  of  as  two  parts:  the  amplitude  [AMP]  on  the  right  hand  side 

(!6)  AMP  =  ha  — 

r 

where  r*  is  the  appropriate  denominator  in  equations  (13,  14,  and  15),  and  the  sinusoidal 
source  phase  O0, 

(17)  Oa=(2  7ft). 

The  difference  in  phase  O,  between  two  locations  is  expressed  by  the  term 

(18)  <”  '  =  -pKr  =  d 

which  is  equal  to  the  exponential  decay  term  (d)  in  equations  (12,  13,  14,  and  15).  Both 
the  amplitude  decay  and  the  degree  of  phase  shift  depend  on  the  ratio  of  hydraulic 
conductivity  to  specific  storage,  which  is  the  hydraulic  diffusivity  (D).  Estimates  of  K 
may  be  inferred  from  equation  (18)  to  compare  with  other  methods  if  Ss  is  assumed. 

The  preceding  equations  can  be  used  to  predict  phase  and  amplitude  versus 
distance  for  homogeneous  systems,  where  K  and  Ss  are  constant.  However,  for 
heterogeneous  systems  where  no  analytical  solutions  are  available,  one  must  resort  to 
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numerical  solutions.  We  postulate  that  relatively  simple  formulas  presented  above  can  be 
used  to  analyze  the  data  for  heterogeneous  cases  by  using  a  distance  weighted  average  for 
the  K.  The  premise  is  that  the  following  replacement  in  the  above  equations  might 
work. 


(19) 


The  index  (I)  indicates  the  present  location  of  r;  so,  the  summation  continues  up  to  the 
present  location  of  r  and  terminates  at  that  point. 

As  indicated  above,  one  must  resort  to  numerical  methods  to  calculate  the  phase 
and  amplitude  relations  with  respect  to  distance  for  heterogeneous  cases  where  K  and  Ss 
change  with  distance.  We  have  developed  numerical  models  for  calculating  the 
amplitude  and  phase  in  the  presence  of  heterogeneity  for  Cartesian,  cylindrical,  and 
spherical  coordinate  systems.  This  research  will  in  later  sections  show  that  the  simple 
replacement  proposed  by  equation  (19),  along  with  equations  (13)  through  (15),  can  be 
used  to  simplify  the  inversion  for  K  in  certain  cases. 

Equations  (14)  and  (15)  represent  the  two  experimental  approaches  utilized  in  this 
research.  The  cylindrical  radial  equation  (14)  describes  the  behavior  of  the  excitation  of 
a  relatively  long  and  small  radius  section  of  screen  that  behaves  as  a  line  source.  Fully 
penetrating  wells  are  often  constructed  at  GEMS.  Any  test  where  the  total  screen  length 
is  excited  is  termed  a  whole  well  test.  The  spherical  radial  equation  (15)  is  a 
representation  of  the  point  source  geometry,  where  the  excited  length  of  well  screen  is 
relatively  short.  To  achieve  this,  either  a  partially  penetrating  well  with  a  relatively  short 
screen  length  or  a  straddle  packer  apparatus  must  be  used.  A  straddle  packer  is  a  double 
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inflatable  packer  arrangement,  which  isolates  a  centralized  interval.  It  is  advantageous  if 
the  packer  apparatus  can  be  deployed  down  typical  2  inch  (5.08  cm)  observation  wells; 
so,  considerable  effort  has  been  expended  to  design  such  packers  for  this  research. 

Previous  studies  have  shown  that  a  line  source  allows  for  higher  energy  input, 
higher  amplitudes,  and  increased  signal  propagation  (Black  and  Kipp,  1981).  A  line 
source  can  create  multiple  ray  paths  to  the  receiver,  decreasing  the  resolution  and  only 
approximating  gross  K  distributions.  High  K  material  can  also  preferentially  propagate 
excess  pore  pressures  generated  by  a  line  source,  which  will  induce  a  vertical  gradient 
and  cross-flow  within  the  aquifer.  Depending  on  the  3-D  heterogeneity  distribution,  this 
cross-flow  will  alter  the  receiver  signal,  similar  to  a  weighted  average,  again  decreasing 
the  resolution.  Even  high  amplitude  line  source  signals  decay  rapidly  in  the  subsurface. 
Most  of  the  decay  is  due  to  the  exponential  term  in  equations  (14)  and  (15).  In  addition, 
the  radial  distance  between  source  and  receiver  wells  will  cause  further  decay  (the 
cylindrical  or  line  source  will  additionally  decay  by  the  inverse  square  root  of  r  [equation 
(14)]  and  the  spherical  or  point  source  will  decay  by  the  inverse  of  r  [equation  (15)]). 
These  additional  amplitude  decay  effects  are  due  to  wavefront  spreading  loss.  However, 
the  point  source  arrangement  may  increase  the  resolution  of  the  K  distribution  profile 
because  of  fewer  ray  path  possibilities. 


The  common  component  of  the  amplitude  decay  and  the  phase  shift  is 


therefore,  it  is  possible  to  compare  the  phase  data  to  the  amplitude  data  (after  correcting 
for  spreading  loss).  Using  aforementioned  assumptions,  estimates  of  K  can  be  obtained 
through  algebraic  manipulation.  However,  this  method  does  not  give  a  specific  value  for 
K,  but  rather  an  average  ratio  of  Ss/K  for  the  signal  travel  path  from  source  well  to 
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receiver  well.  Simple  theory  presented  here  indicates  that  the  phase  and  the  corrected 


amplitude  ratio  should  vary  linearly  with 


and  distance  (r)  from  the  source  well. 


Therefore,  average  parameters  between  well  pairs  may  be  estimated.  Further,  if  multiple 
source  and  receiver  offsets  (relative  to  their  elevations)  are  used,  multiple  diagonal  ray 
paths  may  be  recorded  (Multiple  Offset  Gathers,  MOGs).  This  type  of  testing  is  called 
hydraulic  tomography  (Yeh  and  Liu,  2000;  Bohling  et  al.,  2002;  Bohling  et  al.,  2003), 
and  can  give  more  detailed  information  about  hydraulic  properties  between  wells.  In  the 
first  phase  of  this  project  we  concentrated  on  horizontal  rays  where  the  source  and 
receiver  are  at  the  same  elevation  (Zero  Offset  Profiles,  ZOP).  A  ZOP  survey  is  the 
simplest  tomographical  survey  to  conduct  and  process,  but  can  only  give  information  on 
average  horizontal  aquifer  parameters.  During  follow  up  phases  of  this  project  we  started 
collecting  diagonal  ray  path  data  (MOGs).  These  data  show  the  effects  of  heterogeneity 
in  K  and  offer  the  best  opportunity  to  measure  the  hydraulic  conductivity  distribution. 
Therefore,  we  expended  considerable  effort  trying  to  find  the  optimum  method  of 
processing  these  field  data. 


Field  Methodology 

Recent  studies  at  GEMS  have  utilized  custom-built  straddle  packers  (McElwee 
and  Butler,  1995;  Zemansky  and  McElwee,  2005;  Ross  and  McElwee,  2007),  and 
pneumatic  slug  testing  technique  techniques  (McElwee  and  Zemansky,  2005;  Sellwood, 
2001;  Ross  and  McElwee,  2007).  In  this  work  custom  made  packers  are  used  to  isolate  a 
zone  for  testing.  This  testing  may  either  be  high  resolution  slug  testing  (HRST)  or  cross¬ 
hole  measurement  of  relative  amplitudes  and  phases  for  hydraulic  tomography  . 
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HRST  Techniques 

The  aquifer  material  at  GEMS  exhibits  linear  and  non-linear  responses  to  slug 
testing  (Figure  4).  The  response  of  the  aquifer  material  to  the  slug  can  be  dampened  such 
that  water  levels  in  a  well  return  to  static  head  conditions  with  time  in  a  smooth  non- 
oscillatory  curve.  However,  the  aquifer  can  be  underdamped  and  water  levels  will 
oscillate,  decaying  with  time,  until  pre-test  conditions  are  reached  (Van  Der  Kamp, 

1976).  Theoretical  advances,  presented  by  McElwee  and  Zenner  (1998)  and  McElwee 
(2001,  2002),  have  made  analysis  of  nonlinear  behavior  practical  and  meaningful.  The 
aforementioned  slug  tests  are  localized  tests;  but,  continuous  layers  of  geologic  material 
between  tested  well  pairs  should  correlate  with  HRST  data  from  each  well  in  the  well 
pair. 

Slug  testing  of  an  aquifer  is  an  important  tool  for  determining  aquifer 
heterogeneity  near  a  well.  This  type  of  test  will  average  the  hydraulic  properties  over  a 
limited  volume  of  aquifer.  The  volume  of  aquifer  tested  depends  on  the  length  of  screen 
in  the  aquifer  at  the  tested  well.  A  vertical  profde  of  hydraulic  conductivity  distributions 
can  be  determined  using  high-resolution  slug  testing  in  wells  (Zemansky  and  McElwee, 
2005;  Ross  and  McElwee,  2007),  or  even  with  small  diameter  direct  push  equipment 
(Sellwood,  2001;  Butler  et  al.,  2002a, b;  McCall  et  ah,  2000).  High-resolution  slug  testing 
enables  hydrogeologists  to  examine  vertical  variations  in  K  at  a  much  finer  scale  relative 
to  whole  well  slug  testing. 
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Comparison  of  Tests 
19.5m  Depth 


Comparison  of  Tests 
17.9m  Depth 


Comparison  of  Tests 
15.1m  Depth 


Figure  4.  Three  examples  of  slug  tests  performed  at  GEMS.  Graph  A  displays  no  head 
dependence  and  behaves  linearly.  Graph  B  shows  a  dependence  on  the  initial  slug  height 
and  direction.  Graph  C  is  oscillatory  and  has  some  nonlinear  characteristics. 


The  preferred  method  for  initiating  a  multi-level  slug  test  is  to  use  pneumatics 
(Prosser,  1981;  Zurbuchen  et  al.,  2002;  Zemansky  and  McElwee,  2005;  Ross  and 
McElwee,  2007).  The  advantages  of  using  pneumatics  are  that  nothing  is  added  to  or 
produced  from  the  aquifer  and  less  equipment  is  needed,  which  is  best  for  contaminated 
sites.  The  program  NLSLUG  (McElwee,  2000,  2001)  based  on  the  model  presented  by 
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McElwee  and  Zenner  (1998),  was  used  to  aid  in  the  interpretation  of  oscillatory  and  non- 
oscillatory  hydraulic  head  responses  obtained  from  slug  testing. 

CPT  Techniques 

The  Continuous  Pulse  Test  (CPT)  is  an  exploratory  method  for  extending  slug  test 
results  between  well  pairs  by  propagating  a  sinusoidal  signal.  As  mentioned  earlier,  two 
different  methods  were  used  to  produce  a  sinusoidal  signal:  pneumatic  means  and 
mechanical  pumping.  The  distance  between  wells  in  pairs  tested  and  analyzed  with  the 
CPT  method  in  this  research  have  ranged  from  3  to  1 1 .5  m.  The  instrumentation’s 
ability  to  discern  signal  from  noise  may  be  a  limiting  factor  at  greater  distances.  As  with 
most  geophysical  techniques,  the  equipment  set  up  time  can  consume  considerable  time 
in  the  field.  The  CPT  vertical  profile  method  usually  takes  longer  to  perform  (depending 
on  the  interval  spacing)  than  the  typical  high  resolution  slug  test  vertical  profile  for  a 
given  well. 

In  the  pneumatic  method  an  air  compressor  is  used  to  supply  the  driving  force 
behind  the  CPT  method  and  it  is  connected  to  an  apparatus  attached  to  the  top  of  the 
casing  at  the  well  (Figure  5). 
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Figure  5  a.  The  pneumatic  CPT  equipment  set  up  for  a  line  source  configuration.  A 
signal  generator  opens  and  closes  valves  (Vi  and  V2)  to  control  the  flow  of  air  supplied 
by  the  air  compressor.  The  pressure  transducers  record  the  amplitude  and  phase  at  depth 
Pz  and  a  reference  location  Ps.  This  setup  can  be  easily  modified  for  a  point  source 
configuration  by  using  a  double  packer  to  isolate  the  stressed  interval. 

A  signal  generator  or  computer  controlled  signal  is  used  to  power  servo-controlled  valves 

on  the  apparatus,  which  allows  air  pressure  to  be  increased  in  the  well  or  to  be  released  to 

the  atmosphere.  Increasing  pressure  depresses  the  water  column,  releasing  the  air 

pressure  allows  the  water  column  to  rebound.  A  single  pulse  of  pressure  is  a  slug  test, 

while  stacking  them  one  after  another,  will  create  a  CPT.  The  frequency  and  amplitude 

of  the  CPT  data  should  be  adjusted  to  give  optimal  results  (Engard  et  al.,  2005;  Engard, 

2006).  Figure  5b  shows  the  pneumatic  pressure  manifold  with  the  servo-controlled 

values,  which  was  used  at  the  top  of  the  source  wells. 
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Figure  5b.  Pneumatic  pressure  manifold  with  servo-controlled  valves. 

In  the  mechanical  pumping  method  the  air  is  replaced  by  water  from  a  surface 
reservoir  and  supplied  to  the  aquifer  by  a  pump  through  a  computer  controlled  valve, 
operated  in  such  a  way  that  the  pressure  response  at  the  injection  location  approximates  a 
sine  wave.  This  field  setup  is  shown  in  Figure  5c.  In  this  setup  the  upper  pressure 
transducer  shown  in  Figure  5a  is  not  used.  The  net  result  is  again  a  pressure  signal 
injected  into  the  aquifer  and  measured  by  the  lower  pressure  transducer.  Since  we  are 
continually  injecting  water  in  this  method,  there  is  a  trend  of  increasing  pressure  that 
must  be  removed  by  data  processing  before  the  phase  and  amplitude  are  determined.  In 
an  ideal  setup  the  period  of  this  mechanical  pumping  period  would  be  variable. 

However,  in  this  project  we  were  only  able  to  find  hardware  off  the  shelf  that  would 
allow  a  pumping  period  of  30  seconds.  In  future  research,  it  would  be  preferable  to 
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design  and  build  a  system  that  would  allow  smaller  pumping  periods,  for  reasons  that  will 
be  discussed  later. 


Figure  5c.  Pumped  hydraulic  injection  apparatus  for  the  CPT. 


Surveys  were  done  in  the  form  of  zero  offset  profiles  (ZOP)  and  multiple  offset 
gathers  (MOG).  For  a  ZOP  the  packed-off  source  excitation  interval  with  a  transducer 
and  the  packed-off  receiver  interval  with  a  transducers  are  kept  at  the  same  level,  as  they 
are  moved  through  the  common  screened  interval  of  the  source  and  receiver  wells.  For  a 
MOG,  a  packed-off  source  excitation  interval  with  a  transducer  is  kept  at  a  fixed  depth  in 
the  source  well  while  another  packed-off  receiver  interval  with  a  transducer  is  moved 
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throughout  the  screened  interval  of  the  receiver  well.  For  this  study,  measurements  were 
usually  taken  in  0.30  m  (one  ft)  intervals  (sometimes  1.0  m  or  3  ft  intervals  were  used). 
After  measurements  were  collected  between  one  source  location  and  all  the  receiver 
locations,  the  source  was  moved  by  0.30  m  and  measurements  were  again  collected  at  all 
the  receiver  locations.  The  process  was  repeated  until  rays  had  traveled  from  every 
location  in  the  source  well  to  every  location  in  the  receiver  well  (Figure  6).  The 
collective  examination  of  these  multiple  ray  paths  forms  the  tomographic  study. 

Initially,  a  single-channel  receiver  was  used  in  data  collection.  However,  a  multi¬ 
level  receiver  with  five  pressure  transducers  was  later  constructed  to  expedite  data 
collection.  Pressure  ports  were  located  approximately  1  m  apart  isolated  on  either  side  by 
packers  measuring  approximately  0.6  m  in  length.  The  main  advantage  of  this  apparatus 
is  that  it  allows  efficient  collection  of  multiple  MOGs,  which  are  needed  for  tomographic 
surveys. 

The  MOG  data  taken  from  a  well  pair  should  produce  a  parabolic  phase  shift 
curve  due  to  the  path  lengths  of  the  rays.  Path  lengths  are  greater  for  more  distant  offsets 
(Figure  6).  Larger  phase  and  amplitude  changes  occur  at  these  larger  offsets.  If  the 
source  is  in  the  middle  of  the  well,  the  greatest  distance  and  therefore  greatest  change  in 
amplitude  and  phase  should  occur  when  the  receiver  is  at  the  top  or  bottom.  The  shortest 
distance  is  when  the  source  and  receiver  are  at  the  same  depth.  The  general  shape  should 
be  a  parabola  with  distortions  due  to  heterogeneity.  When  the  source  is  at  the  top,  the 
shortest  distance  is  to  the  receiver  location  at  the  same  depth  and  the  greatest  distance  is 
to  the  receiver  location  at  the  bottom  of  the  well.  The  curve  should  therefore  have  a  half¬ 
parabola  shape  when  the  source  is  at  the  top  of  the  well.  The  same  is  true  when  the 
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source  is  at  the  bottom  of  the  well.  Examples  of  these  parabolic  shapes  are  shown  in 


Figures  6. 


Sou  rce  we  1 1  Receiver  we  1 1 


Figure  6.  MOG  setup  for  the  tomographic  study. 

Pressure  transducers  were  used  to  monitor  pressure  head  fluctuations  in  both  the 
source  well  and  at  the  observation  wells.  The  data  were  collected  from  the  pressure 
transducers  by  a  data-logger  and  stored  on  a  field  computer  for  later  analysis.  Data  were 
typically  recorded  at  a  20  Hz  sampling  rate,  which  provided  sufficient  temporal 
resolution.  The  field  computer  and  data  logger  allowed  real-time  monitoring  of  the  CPT 
records. 
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Vertical  Sensor  Array 

Throughout  the  project  we  have  continued  to  improve  the  design  of  the  vertical 
sensor  array.  Moving  the  receiver  location  to  many  discrete  locations  along  the  receiver 
well  screen  is  very  time  consuming.  To  speed  this  process,  we  designed  a  vertical  sensor 
array  with  five  pressure  transducers  and  six  packers.  Each  transducer  is  isolated  by 
packers  above  and  below,  to  allow  measurements  to  be  made  on  a  0.3  m  (1  ft)  section  of 
the  receiver  well  screen.  The  transducers  are  located  every  0.91  m  (3  ft)  along  the  array, 
with  0.6  m  (2  ft)  length  packers  between.  The  array  may  be  moved  up  in  0.3  m  (1  ft) 
increments  two  times  to  allow  uniform  coverage  of  the  first  section  of  the  screen  at  0.3  m 
(1  ft)  increments.  Nearly  complete  coverage  of  the  1 1  m  screen  can  be  achieved  by 
pulling  the  vertical  sensor  array  3.9  m  (13  ft)  and  repeating  the  sequence  described  above. 
In  this  way  recording  six  records  with  the  vertical  sensor  array  is  equivalent  to  30  records 
with  the  single  receiver  setup.  This  increases  the  speed  of  data  collection.  Pictures  of  the 
first  generation  vertical  sensor  array  are  shown  below  in  Figure  7a. 

According  to  the  project  plan,  we  were  to  adapt  the  sinusoidal  source  to  be  used 
with  a  Geoprobe  unit.  Initially,  the  Geoprobe  unit  is  used  to  advance  a  drill  string  to  the 
bottom  of  the  aquifer.  At  that  point,  the  drive  tip  is  replaced  with  a  source  unit  capable  of 
generating  a  sinusoidal  signal  and  then  retracted  in  stages  to  occupy  each  desired  source 
location.  During  data  collection  with  the  source  on  the  end  of  the  Geoprobe  unit  we  used 
two  vertical  sensor  arrays,  since  it  was  expensive  to  have  the  Geoprobe  on  site.  We  took 
this  opportunity  to  improve  the  design  of  the  vertical  sensor  array  and  to  construct 
another  one.  This  improved  design  is  shown  in  Figure  7b.  The  basic  dimensions  are  the 
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same;  however,  the  port  design  and  method  of  connection  of  the  packers  was  changed  to 
allow  easier  deployment  and  retrieval. 


Figure  7a.  First  generation  vertical  sensor  array. 
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New  Wells  Installed 

In  late  October  2007,  three  wells  were  added  to  GEMS.  The  wells  were  chosen  to 
provide  better  coverage  of  the  area  under  study  by  hydraulic  tomography.  The  wells 
were  installed  using  the  direct  push  method  with  a  Geoprobe  unit  from  the  Kansas 
Geological  Survey.  The  wells  initially  installed  for  this  project  were  HT-1,  HT-2,  and 
HT-3.  The  new  wells  are  HT-4,  HT-5,  and  HT-6.  All  of  these  wells  and  others 
previously  used  for  hydraulic  tomography  work  are  shown  below  in  Figure  8.  After 
installation  and  development,  the  wells  were  surveyed  to  establish  the  elevation  of  the  top 
of  each  casing.  The  Geoprobe  source  well  location  is  also  shown.  Various  radii  between 
wells  were  measured  for  future  analysis  of  the  cross- well  data.  All  of  this  information 
about  the  various  wells  that  were  candidates  for  tomographic  study  is  shown  in  Table  1. 
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Figure  8.  Relative  well  locations  at  GEMS  (north  is  up).  This  shows  the  locations  of  the 
new  wells  installed  in  Oct.  2007  and  the  Geoprobe  source  location,  in  addition  to  older 
wells  previously  used  in  this  study. 
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Table  1.  Well  Information 


Location 

Elevation  ft 

Elevation  m 

Depth  ft 

Depth  m 

Screen  ft 

Screen  m 

827.556 

252.239 

830.005 

252.986 

72.3 

10.67 

829.66 

252.880 

72.4 

10.67 

■sa 

829.705 

252.894 

~70. 

10.67 

830.129 

253.023 

72.2 

22.01 

10.67 

829.651 

252.878 

71.9 

21.92 

10.67 

830.272 

253.067 

~72. 

-21.9 

10.67 

7-1 

828.342 

252.479 

68.85 

20.99 

30.0 

9.14 

11-1 

828.358 

252.484 

69.40 

21.16 

45.0 

13.72 

Inj.  Well 

829.794 

252.921 

71.09 

21.67 

34.0 

10.36 

HT-GP  ref. 

828.82 

252.62 

Well  to  Well  Radial  Distances,  r 


Well 

Well 

Radius  (m) 

Radius  (ft) 

HT-3  to 

HT-1 

4.77 

15.65 

HT-3  to 

HT-2 

4.36 

14.31 

HT-3  to 

HT-4 

4.46 

14.62 

HT-3  to 

HT-5 

4.21 

13.81 

HT-3  to 

HT-6 

3.99 

13.10 

HT-3  to 

HT-GP 

4.25 

13.94 

HT-2  to 

HT-GP 

4.23 

13.88 

HT-6  to 

7-1 

2.70 

8.85 

HT-6  to 

11-1 

7.19 

23.58 

HT-6  to 

Inj.  Well 

4.04 

13.26 

Inj.  Well  to 

HT-1 

4.28 

14.05 

Inj.  Well  to 

HT-4 

8.67 

28.45 

Inj.  Well  to 

HT-5 

11.55 

37.89 

Inj.  Well  to 

HT-2 

11.49 

37.70 

Inj.  Well  to 

HT-3 

7.66 

25.15 

7-1  to 

HT-2 

6.94 

22.79 

7-1  to 

HT-5 

9.18 

30.10 

7-1  to 

HT-3 

5.13 

16.84 

7-1  to 

HT-4 

9.00 

29.53 

7-1  to 

HT-1 

6.46 

21.20 

HT-6  to 

HT-1 

3.79 

12.42 

HT-1  to 

HT-4 

4.40 

14.44 

HT-4  to 

HT-5 

4.63 

15.21 

HT-5  to 

HT-2 

4.57 

15.00 

HT-2  to 

HT-6 

7.40 

24.28 

32 


Data  Processing,  Modeling,  and  Inversion 

Data  Processing 

Data  processing  for  the  ZOP  data  was  done  with  FitAmpPhaseV8,  a  program 
written  in  Visual  Basic  by  Carl  McElwee.  The  program  fits  sine  waves  to  the  transducer 
data  taken  in  the  field  and  generates  plots  of  the  amplitude  ratio  and  phase  shift  (x-axis) 
between  the  source  and  receiver  transducers.  All  values  are  plotted  against  location  (y- 
axis).  The  program  analyzes  data  for  a  single  source  location  at  a  time.  For  each  MOG, 
the  amplitude  ratio  and  phase  shift  between  the  two  (upper  and  lower)  source  transducers 
should  plot  as  a  vertical  line,  as  there  is  no  change  in  material  within  the  source  well 
itself.  The  amplitude  ratio  and  phase  shift  between  the  (lower)  source  and  receiver 
transducers  should  both  plot  as  nearly  parabolas  or  half-parabolas.  If  the  source  location 
is  near  the  middle  of  the  well,  the  shape  will  be  a  full  parabola,  and  the  shape  will  only  be 
half  a  parabola  if  the  source  is  near  either  the  top  or  bottom  of  the  well.  The  shape  should 
be  nearly  parabolic  assuming  no  change  in  aquifer  material  (or  measurement  error),  so 
any  deviations  from  the  overall  parabola  must  be  due  to  changes  in  K  (or  some 
experimental  error). 

Data  processing  for  the  MOG  data  was  done  with  FitAmpPhaseVlOHT  (early  3-4 
sec.  data)  or  FitAmpPhaseV12HT  (later  30  sec.  and  3  sec.  data).  These  versions  of  the 
program  analyze  all  five  receiver  transducers  at  once.  Aside  from  accounting  for 
multiple  receiver  transducers,  the  programs  are  based  on  the  same  algorithms  as 
FitAmpPhaseV8.  Some  improvements  or  changes  were  made  along  the  way,  giving  rise 
to  version  numbers.  Version  10  had  to  fit  the  period  as  a  parameter,  because  the  early 
data  used  a  frequency  generator  for  manual  period  selection  (which  was  never  quite  the 
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same).  Version  12  deals  with  data  whose  period  is  computer  controlled  and  consistent 
from  record  to  record,  and  also  removes  any  trend  of  increasing  pressure  due  to  the 
continual  injection  for  the  mechanical  pumping  system.  Plots  are  generated  for  the 
receiver  location  versus  both  amplitude  and  phase  shift.  The  raw  data  and  the  fitted  sine 
wave  for  a  single  receiver  location  are  shown  below  in  Figure  9  for  some  example  data. 
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Figure  9:  The  data  for  one  particular  receiver  location  in  the  FitAmpPhaseVlOHT 
program.  Three  plots  are  shown:  one  plot  for  each  of  the  two  source  transducers  and  one 
plot  for  a  receiver  transducer.  The  raw  data  are  shown  in  blue  while  the  fitted  sine  wave 
is  shown  in  pink. 
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High  resolution  slug  test  (HRST)  data  were  processed  using  the  program 
NLSLUG  (McElwee,  2000),  developed  by  Carl  McElwee  using  Fortran  and  run  from 
Microsoft  Excel.  Water  and  air  pressure  transducers  are  used  to  record  the  initial  height 
of  the  slug  test.  For  each  record,  a  time  break  is  chosen  to  begin  measuring  time,  and 
static  values  at  long  times  are  determined  for  a  base  line.  Multiple  initial  heads  are  used. 
If  the  results  are  independent  of  initial  head  and  behave  linearly,  all  records  lie  on  top  of 
each  other.  Usually  the  records  do  not  completely  overlie  one  another,  so  there  can  be 
problems  with  both  directionality,  i.e.  positive  or  negative  initial  head,  and  repeatability. 
Mobile  fine  sediments  could  explain  both  problems  (McElwee  and  Zemansky,  2005). 
Slug  testing  can  cause  fine  sediments  to  move,  and  these  sediments  may  move  more 
easily  into  the  well  than  away  from  the  well,  creating  an  annulus  containing  more  fine 
material  at  some  radius.  HRST  data  for  the  wells  in  this  study  were  processed  by  Brett 
Engard  and  Pema  Deki  (Appendix  B).  The  HRST  results  can  be  used  to  constrain  the 
inversion  to  ensure  that  the  interwell  K  values  remain  in  the  range  observed  in  HRST 
results. 

Straight  Ray  Modeling  With  Spatially  Weighted  K  Values 

Typical  hydraulic  tomography  inversions  use  nonlinear  least  squares  fitting,  a 
numerical  model,  and  iterations  to  get  the  best  fit,  a  process  that  can  take  much  time  and 
computing  power.  An  approximation  for  the  numerical  model  has  been  used  in  this 
research  using  straight,  spatially  weighted  ray  paths.  The  path  length  in  each  zone  of 
differing  K  is  multiplied  by  a  coefficient  involving  K  to  get  the  phase.  This  is  a  direct 
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implementation  of  equation  (19),  and,  the  accuracy  of  this  approximation  must  be 
checked  by  comparison  to  a  numerical  model  (which  will  be  done  in  a  later  section). 

Ray  path  data  were  generated  by  HydraulicTomAnal(V19  and  V21,  slight 
differences  in  version  numbers),  developed  by  Carl  McElwee  in  Microsoft  Excel.  The 
field  area  was  divided  into  a  grid  system  with  approximately  evenly  spaced  divisions  in 
the  horizontal  and  vertical  directions.  Each  box  within  the  grid  is  referred  to  as  an 
element.  The  model  was  divided  into  a  series  of  nodes,  elements,  and  grid  spaces  (Figure 


Figure  10:  Depiction  of  a  node,  an  element,  and  a  grid  space. 


Nodes  are  any  of  the  individual  points  throughout  the  grid.  The  vertical  or  horizontal 
spaces  between  two  nodes,  Ax  and  Az,  are  known  as  grid  spacings.  An  element  is  the 
rectangular  area  enclosed  by  four  adjacent  node  points. 
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The  values  of  K  may  be  given  either  by  nodes  or  elements.  If  K  is  given  by  nodes 
the  K  value  is  linearly  interpolated  between  nodes.  If  K  is  given  by  elements  then  K  is 
constant  everywhere  within  that  element.  The  program  allows  zones  of  constant  K  by 
grouping  elements  or  nodes  together  (depending  on  how  K  is  specified),  which  share  a 
common  value  of  K.  The  program  computes  the  distance  of  each  ray  path  through  every 
element  based  on  the  Pythagorean  Theorem.  The  phase  shift  for  each  element  is  given  by 
an  implementation  of  equation  (19),  multiplying  the  path  length  by  the  appropriate  value 
of  K.  Path  lengths  through  each  element  or  zone  and  phase  shift  values  for  each  total  ray 
path  are  then  available  as  output  for  use  in  an  inversion  program. 

Data  Inversion  for  K  Values 

Once  we  have  obtained  data,  either  from  the  field  or  by  running  a  numerical 
model,  we  have  a  time  series  representing  the  source  at  some  location  and  a  receiver  at 
another  location.  These  time  series  may  be  analyzed  to  find  the  phase  and  relative 
amplitude  at  the  receiver  by  the  method  outlined  above.  We  may  have  a  few  (ZOP)  or 
many  (MOG)  ray  paths  to  analyze.  In  the  case  of  synthetic  numerical  data  we  know  the 
input  K  values  and  can  test  the  inversion  process  to  see  how  well  these  values  are 
returned.  In  the  case  of  real  field  data  we  do  not  know  the  correct  values.  In  any  case, 
we  must  assume  some  model  structure  for  the  data  in  order  to  perform  the  inverse.  This 
will  involve  assuming  the  number  of  nodes,  elements,  and  zones  of  constant  K  to  use.  It 
is  well  known  in  the  inverse  literature  that  the  inversion  may  be  non-unique  due  to  a 
number  of  factors,  including  model  structure  and  measurement  error.  In  a  later  section 
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the  stability  and  uniqueness  of  various  models  will  be  investigated.  In  this  section  we 
merely  describe  the  inversion  program  that  will  be  used. 

The  output  of  the  program  HydraulicTomAnal(V19  or  V21)  contains  path  lengths 
through  each  K  zone  and  total  phase  shift  for  each  straight  ray  in  the  data  set.  These  data 
may  be  transferred  to  the  LeastSquareSVDV13  or  VI 5  (slight  program  differences  in 
version  numbers)  program,  developed  by  Carl  McElwee  in  Microsoft  Excel.  The  SVD, 
or  Singular  Value  Decomposition,  program  performs  a  least  squares  fitting  inversion 
using  zone  path  lengths  and  total  phase  shift  values  for  all  ray  paths,  to  obtain  K  values 
for  each  zone  by  using  a  set  of  linear  equations  (Aster  et  al.,  2005).  Equations  used  in  the 
program  do  not  require  iterations  with  a  numerical  model  because  they  are  linear  due  to 
the  straight  ray  approximation  used  here.  The  SVD  method  divides  G  (matrix  of  zone  ray 
path  lengths),  an  m  (number  of  ray  paths  and  equations)  by  n  (number  of  zones  and 
unknowns)  matrix  into  the  following  equation: 

G  =  UWVT  (20) 

where  U  is  an  m  by  m  orthogonal  matrix,  W  is  an  m  by  n  matrix  with  nonnegative 
diagonal  elements  known  as  singular  values,  V  is  an  n  by  n  orthogonal  matrix,  and  the  T 
superscript  indicates  that  V  is  a  transpose  matrix.  Standard  deviations  of  the  fitted 
parameters  (K)  are  calculated  based  on  the  goodness  of  fit.  The  program  has  the  added 
feature  of  reducing  to  deterministic  inversion  if  the  number  of  equations  (ray  paths)  and 
unknowns  (K  values)  are  equal,  provided  the  matrices  are  non-singular.  This  is  the 
program  we  will  use  to  theoretically  evaluate  model  stability  and  uniqueness  and  to 
perform  inversion  of  the  field  data  to  obtain  K  distributions. 
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The  inversion  program  (LeastSquareSVDV13  or  VI 5)  has  options  which  can  be 
selected.  If  synthetic  data  are  being  inversed,  then  the  program  can  be  run  in  Monte 
Carlo  mode  with  a  specified  per  cent  of  random  error  added  to  the  data,  and,  a  large 
number  of  inversions  may  be  analyzed  statistically.  The  program  allows  an  initial 
estimate  of  the  K  values  to  be  made  and  a  penalty  constraint  to  be  applied  as  the  inferred 
values  deviate  from  the  initial  estimate.  This  is  a  valuable  tool  to  constrain  any  non¬ 
uniqueness  tendencies,  if  one  has  independent  data  regarding  K.  In  our  case  we  have 
some  K  values  from  HRST  which  can  be  used  as  constraints. 

Finite  Difference  Numerical  Modeling 

As  indicated  above,  one  must  resort  to  numerical  methods  to  calculate  the  phase 
and  amplitude  relations  with  respect  to  distance  for  heterogeneous  cases  where  K  and  Ss 
are  changing  with  distance.  We  have  adapted  a  numerical  model  written  by  Carl 
McElwee  during  his  years  of  teaching  groundwater  modeling  at  the  University  of  Kansas. 
This  model  allows  calculation  of  the  pressure  response  at  an  arbitrary  receiver  location  in 
response  to  an  oscillatory  input  at  any  given  finite  location  on  the  well  screen  in  the 
presence  of  heterogeneity  for  Cartesian  and  radial  coordinate  systems. 

GSIT2DTVHeter  is  a  computer  program  in  VBA  (Visual  Basic  for  Applications) 
written  within  the  Excel  spreadsheet  environment  to  evaluate  the  Cartesian  or  radial 
equation  for  a  two-dimensional  system  of  coordinates. 

ah  Q 

V  •  {KWh)  =  Ss  —  -  ,  0  <  x  or  r  <  MaxX,  0  <  z  <  MaxZ 

The  above  equation  has  been  written  for  groundwater  flow.  Hydraulic  head  (Length)  is 

3 

given  by  h,  the  sources  or  sinks  of  water  (Q)  for  the  system  (Length  /Time)  represents 
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water  added  to  (positive)  or  subtracted  from  (negative)  the  system  per  unit  volume  (AV), 
K  is  the  hydraulic  conductivity  (Length/Time),  and  Ss  is  the  specific  storage  (1/Length). 
The  physical  parameters  Kx  (hydraulic  conductivity  in  x  direction),  Kz  (hydraulic 
conductivity  in  z  direction),  Ss  (specific  storage),  and  m  (aquifer  thickness)  may  be 
specified  for  every  node,  simulating  heterogeneity.  MaxX  and  MaxZ  are  the  lengths  in 
the  x  and  z  directions  of  the  region  for  solution.  Boundary  conditions  on  the  four  sides 
can  be  either  head  specified  or  derivative  specified.  For  the  radial  case  an  alternate  is  to 
specify  the  pump  rate  of  the  well  at  the  left-hand  boundary  (assumed  well  screen 
location);  in  this  case,  the  boundary  condition  on  the  well  screen  can  vary  from  node 
point  to  node  point.  This  feature  allows  us  to  put  a  sinusoidally  varying  pumping  rate  at 
one  or  more  nodes  and  then  specify  a  barrier  boundary  condition  (simulating  the  presence 
of  a  packer)  for  the  other  well  screen  nodes. 

The  use  of  GSIT2DTYHeter  allows  us  to  generate  synthetic  field  records  for 
various  source  and  receiver  locations,  simulating  ZOP  and  MOG  field  records.  These 
simulated  field  records  may  be  processed  thorough  the  FitAmpPhase  program  to  obtain 
the  relative  amplitude  and  phase  shift  at  the  receiver  and  then  processed  further  by  the 
Fly draulicTom Anal  and  LeastSquareSVDY  programs  as  discussed  above  in  order  to  test 
the  validity  of  the  straight,  spatially  weighted  ray  path  approximation. 

Investigation  of  Straight  Ray  Approximation 

As  presented  earlier  in  this  report  it  is  postulated  that  perhaps  a  spatially  weighted 
average  K  value  [equation  (19)]  could  be  substituted  into  the  homogeneous  analytical 
solution  [equations  (13)  through  (15)]  as  an  approximation  to  the  heterogeneous  case. 
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Since  no  analytical  solutions  exist  for  arbitrary  heterogeneous  systems,  we  must  resort  to 
numerical  modeling  to  check  this  approximation.  Modeling  studies  were  performed  to 
compare  results  from  the  spatially  weighted  ray-tracing  method  with  those  from  a 
numerical  model  for  Cartesian  and  cylindrical  coordinates.  The  goal  is  to  see  if  the 
simple  replacement  proposed  above  can  simplify  the  inversion  for  K. 

Using  the  output  of  the  numerical  models,  we  used  an  early  version  of  the 
FitAmpPhase  program  to  calculate  the  phase  and  amplitude  as  a  function  of  distance  for 
heterogeneous  models  for  the  Cartesian  and  cylindrical  coordinate  systems  (no  variation 
in  the  vertical  direction).  Using  equation  (19)  as  an  approximation  we  can  calculate  the 
K  values  expected  from  these  values  of  phase.  We  looked  at  systems  consisting  of  blocks 
of  material  with  differing  K  and  for  systems  where  K  varied  systematically,  such  as  in  a 
linear  trend.  As  can  be  seen  from  the  data  presented  in  Table  2,  the  agreement  between 
the  numerical  data  and  the  theory  using  a  spatially  weighted  average  to  solve  for  K  is 
excellent,  except  near  boundaries  and  near  the  origin.  The  calculated  values  for  K  were 
determined  by  considering  the  phases  from  the  numerical  models.  The  results  for  K 
using  the  amplitude  data  are  similar  but  have  a  little  more  error  near  the  origin.  We 
believe  this  technique  will  work  for  the  spherical  coordinate  system  also  (allowing 
variation  in  the  vertical  direction)  and  is  the  subject  of  a  following  section.  This 
simplification  in  solving  for  K  should  make  the  tomographic  inversion  considerably 
simpler  than  if  a  full  numerical  model  was  needed  to  solve  for  K. 
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Table  2.  Comparison  of  approximate  results  for  hydraulic  conductivity  compared  with 
true  numerical  model  values. 


Cartesian  coordinates: 

Two  zones  for  K 


X 

0 

5 

10 

15 

20 

25 

30 

35 

40 

Amplitude 

1 

0.333887 

0.111552 

0.03766 

0.010765 

0.0046 

0.002118 

0.000974 

0.000449 

Phase 

0 

-0.17316 

-0.34662 

-0.51784 

-0.68759 

-0.81992 

-0.94258 

-1.0656 

-1.18951 

Cal.  K 

0.002985 

0.002968 

0.003198 

0.002403 

0.005955 

0.005943 

0.005891 

0.005786 

True  K 

0.003 

0.003 

0.003 

0.003 

0.006 

0.006 

0.006 

0.006 

Linearly  varying  K 


X 

0 

5 

10 

15 

20 

25 

30 

35 

40 

Amplitude 

1 

0.33608 

0.126081 

0.051371 

0.022336 

0.010236 

0.004895 

0.002428 

0.001256 

Phase 

0 

-0.16353 

-0.3115 

-0.44764 

-0.5744 

-0.69342 

-0.806 

-0.91361 

-1.01703 

Cal.  K 

0.003653 

0.004393 

0.005133 

0.005875 

0.006625 

0.007354 

0.00797 

0.008727 

True  K 

0.0036 

0.00435 

0.0051 

0.00585 

0.0066 

0.00735 

0.0081 

0.00885 

Cylindrical  coordinates: 

Two  zones  for  K 


r 

0.0833 

1.0231 

5.1071 

10.0331 

15.2399 

20.3548 

25.4931 

30.9181 

35.1619 

39.9883 

Amplitude 

1.0000 

0.3834 

0.0805 

0.0200 

0.0053 

0.0013 

0.0005 

0.0002 

0.0001 

0.0000 

Phase 

0.0000 

-0.0494 

-0.2012 

-0.3753 

-0.5565 

-0.7248 

-0.8690 

-1.0028 

-1.1080 

-1.2289 

Cal.  K 

0.0013 

0.0026 

0.0029 

0.0030 

0.0033 

0.0045 

0.0059 

0.0058 

0.0057 

True  K 

0.0030 

0.0030 

0.0030 

0.0030 

0.0030 

0.0060 

0.0060 

0.0060 

0.0060 

Linearly  varying  K 


r 

0.0833 

1.0231 

5.1071 

10.0331 

15.2399 

20.3548 

25.4931 

30.9181 

35.1619 

39.9883 

Amplitude 

1.0000 

0.3696 

0.0766 

0.0211 

0.0068 

0.0025 

0.0010 

0.0004 

0.0002 

0.0001 

Phase 

0.0000 

-0.0466 

-0.1861 

-0.3334 

-0.4757 

-0.6052 

-0.7271 

-0.8484 

-0.9394 

-1.0395 

Cal.  K 

0.0026 

0.0035 

0.0044 

0.0052 

0.0060 

0.0068 

0.0074 

0.0081 

True  K 

0.0031 

0.0038 

0.0045 

0.0053 

0.0060 

0.0068 

0.0076 

0.0083 

As  shown  above,  the  homogeneous  equations  can  be  used  to  predict  K  based  on 
the  measurable  amplitude  decay  and  phase  shift.  However,  the  values  obtained  for  the 
horizontal  rays  must  be  interpreted  as  spatially  weighted  averages  over  the  horizontal 
distance  between  wells.  Equations  (14)  and  (15)  represent  the  two  experimental 
approaches  utilized  in  this  research.  The  cylindrical  radial  equation  (14)  describes  the 
behavior  of  the  excitation  of  a  relatively  long  and  small  radius  section  of  screen  and  is 
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considered  to  behave  like  a  line  source.  Fully  penetrating  wells  are  often  constructed  at 
GEMS.  Any  test  where  the  total  screen  length  is  excited  is  termed  a  whole  well  test.  The 
spherical  radial  equation  (15)  is  a  representation  of  the  point  source  geometry,  where  the 
excited  length  of  well  screen  is  relatively  short.  To  achieve  this,  either  a  partially 
penetrating  well  with  a  relatively  short  screen  length  or  a  straddle  packer  apparatus  must 
be  used. 

Finally,  we  investigate  the  validity  of  the  spatially  weighted  straight  ray 
approximation  where  vertical  variation  can  occur  in  K  and  rays  are  allowed  to  be 
diagonal  in  addition  to  horizontal  (spherical  geometry  case).  Again,  modeling  studies 
were  preformed  to  compare  results  from  the  spatially  weighted  ray  tracing  method  with 
those  from  a  numerical  model.  The  numerical  model  and  straight  ray  method  were  both 
used  to  simulate  the  phase  shift  of  108  rays  between  a  theoretical  well  pair  with  three 
CPT  source  locations,  each  with  36  corresponding  receiver  locations.  Modeling  was 
completed  for  both  the  3 -sec  and  30-sec  CPTs  to  compare  the  difference  between  the  two 
source  methods.  The  aquifer  between  the  well  pair  was  simulated  by  a  3-element,  8- 
node,  model  which  corresponds  to  the  screen  interval  [10.68  m  (35  ft)]  and  radial 
distance  [5.85  m  (19.20  ft)]  between  the  theoretical  well  pair.  The  upper,  middle  and 
lower  elements  are,  respectively,  4.88  m  (16  ft),  0.92  m  (3  ft),  and  4.88  m  (16  ft)  thick. 
The  upper,  middle  and  lower  elements  have  K  values  of  0.0009  m/sec  (0.003  ft/sec), 
0.0018  m/sec  (0.006  ft/sec),  and  0.0009  m/sec  (0.003  ft/sec),  respectively  (Fig.  11).  A 
representative  Ss  value  of  0.00018  was  also  assumed  for  the  verification  modeling. 
Although  these  values  were  arbitrarily  chosen,  they  fall  within  the  range  of  values 
observed  at  GEMS  and  are  consistent  with  Wachter’s  (McElwee  et  al.,  2007;  Wachter, 
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2008;  Wachter  et  al.,  2008)  earlier  verification  of  the  heterogeneity  extension  using  a  4- 
sec  pneumatic  CPT.  The  numerical  phase  data  from  this  model  comprise  a  theoretically 
perfect  CPT  data  set  and  phase  data  from  the  straight  ray  model  should  closely 
approximate  it.  However,  the  numerical  model  does  use  a  barrier  boundary  on  the  top 
and  bottom  rather  than  an  infinite  domain,  so  some  boundary  effects  are  expected.  In  any 
case,  good  agreement  between  the  two  methods  is  a  line  of  evidence  supporting  the 
heterogeneous  extension  [equation  (19)]  adapted  for  this  research. 

Wachter’ s  4-sec  CPT  numerical  validation  of  the  heterogeneity  extension  was 
reproduced  with  the  latest  version(s)  of  the  Visual  Basic  data  processing  programs  so  his 
verification  could  be  compared  to  the  numerical  verification  of  the  30-sec  CPT  data  used 
in  the  later  stages  of  this  research.  Numerical  modeling  simulated  three  MOG  data  sets 
from  source  locations  at  0.305-0.610  m  (1-2  ft),  5.486-5.791  m  (18-19  ft),  and  10.668- 
10.973  m  (35-36  ft),  which  correspond  to  the  lower,  middle,  and  upper  intervals  of  the 
aquifer  model  (Fig.  11-13).  The  numerical  model  had  36  rows  to  simulate  each  of  the 
36  theoretical  receiver  locations  in  a  MOG.  To  simulate  a  file  of  head  data  from  the  CPT 
source  and  receiver  transducers,  numerical  phase  data  was  parsed  from  the  numerical 
model  rows  at  radial  distances  which  correspond  to  the  center  of  the  source  and  receiver 
well  locations  (e.g.,  0.25  m  [0.833  ft])  and  5.85  m  [19.20  ft])  and  were  saved  to  a  text 
file.  FitAmpPhase  used  the  text  files  to  calculate  the  numerical  phase  shift  for  each  of  the 
MOGs.  HydraulicTomAnal  was  used  to  create  an  element  matrix  of  the  aquifer  and 
apply  the  straight  ray  approximation  method  through  the  matrix  to  generate  the  straight- 
ray  phase  shift  data  for  all  three  MOGs.  The  element  matrix  was  imported  into 
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LeastSquaresSVD  and  both  numerical  and  straight-ray  phase  shift  data  for  the  three 
MOGs  were  inverted  through  the  element  matrix  to  calculate  diffusivity  and  solve  for  K. 

The  current  version  of  the  SVD  inversion  program  also  has  the  ability  to  perform 
Monte  Carlo  simulations  using  random  error,  rather  than  running  individual  simulations. 
Monte  Carlo  simulations  were  run  with  both  5%  and  10%  random  noise  for  1000 
simulations.  The  5%  random  noise  approximates  the  expected  variation  in  the  field  due 
to  instrument  imprecision  and  ambient  noise  and  the  10%  random  noise  simulates  the 
expected  worse-case  scenario  of  signal  inference.  Verification  of  the  heterogeneity 
extension  and  comparison  of  the  4  and  30-sec  CPT  sources  are  further  discussed  below. 

The  numerical  phase  shift  from  the  4  sec  sources  were  compiled  and  compared  to 
their  corresponding  straight-ray  approximation  to  evaluate  the  relative  goodness  of  fit 
between  the  simulated  field  data  and  its  model  approximation  (Fig.  1 1  -  13).  Because  the 
phase  originates  from  synthetic  data,  the  two  curves  should  fit  relatively  close.  The  4-sec 
CPT  phase  shift  values  from  the  spatially  weighted  ray  method  and  the  numerical  model 
for  the  upper,  middle,  and  lower  source  locations  were  in  good  agreement  with  each  other 
except  for  some  slight  boundary  effects  (Fig.  11-13).  There  was  some  deviation  of  the 
straight  ray  phase  shift  at  the  middle  source  location  through  the  thinner,  middle  layer 
(Fig.  12).  Straight  rays  projected  through  this  element  more  directly  measure  the  K 
without  the  averaging  across  the  middle  layer  from  the  numerical  model  due  to 
wavelength  considerations  and  result  in  the  higher  K  values  (i.e.,  low  phase)  seen  in  this 
layer  of  the  graph.  Overall,  the  data  fit  is  good  indicating  resolution  of  about  1  m  (3  ft) 
layers  with  a  4-sec  period,  reconfirming  Wachter’s  (McElwee  et  al.,  2007;  Wachter, 

2008;  Wachter  et  al.,  2008)  assessment  of  the  resolution. 
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Numerical  Model  Phase  Shift  vs 
Spatially  Weighted  Straight  Ray  Phase  Shift 
4  Second  CPT  Period 

0.0-0.305  m  (0-1  ft)  Theoretical  Source  Location 
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♦  Numerical  Phase  Shift  from  FitAmpPhase  — ■ - Straight  Ray  Approximation  from  HydraulicTomAnal  •  Source  Location 


Figure  1 1  -  A  comparison  of  4  second  CPT  period  phase  shift  values  from  a  numerical 
model  and  the  spatially  weighted  ray  path  method  at  the  0.0-0.305  m  (0-1  ft)  source 
location. 


Numerical  Model  Phase  Shift  vs 
Spatially  Weighted  Straight  Ray  Phase  Shift 
4  Second  CPT  Period 

5.19-5.49  m  (17-18  ft)  Theoretical  Source  Location 


Phase  Shift 

1.6  1.4  1.2  1  0.8  0.6  0.4  0.2  0 


♦  Numerical  Phase  Shift  from  FitAmpPhase  — ■ - Straight  Ray  Approximationfrom  HydraulicTomAnal  •  Source  Location 


Figure  12  -  A  comparison  of  4  second  CPT  period  phase  shift  values  from  a  numerical 
model  and  the  spatially  weighted  ray  path  method  at  the  5.19-5.49  m  (17-18  ft)  source 
location. 
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Numerical  Model  Phase  Shift  vs 
Spatially  Weighted  Straight  Ray  Phase  Shift 
4  Second  CPT  Period 

10.37-10.68  m  (34-35  ft)  Theoretical  Source  Location 


Phase  Shift 
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♦  Numerical  Phase  Shift  from  FitAmp Phase  — ■ - Straight  Ray  Approximation  from  HydraulicTomAnal  •  Source  Location 


Figure  13  -  A  comparison  of  4  second  CPT  period  phase  shift  values  from  a  numerical 
model  and  the  spatially  weighted  ray  path  method  at  the  10.37-10.68  m  (34-35  ft)  source 
location. 

The  30-second  CPT  phase  shift  values  from  the  spatially  weighted  ray  method 
and  the  numerical  model  at  the  upper,  middle,  and  lower  source  locations  were  in 
reasonable  agreement  although  the  data  resolution  or  overlap  of  the  two  curves  was  not  as 
precise  as  the  4-second  MOG  data  sets.  The  resolution  of  a  longer  period  signal  is 
expected  to  be  less  due  to  the  longer  wavelength  of  the  propagating  signal  (therefore 
averaging  over  a  larger  volume)  and  results  such  as  this  are  a  piece  of  evidence  to  support 
that  theory.  In  general,  the  data  curves  are  similar  and  the  slight  boundary  effects  are  still 
present  (Fig.  14  -  16).  Again,  there  was  some  deviation  of  the  straight  ray  phase  shift 
through  the  thinner,  middle  layer  (Fig.  15).  Also,  the  two  phase  shift  curves  were  offset 
slightly  at  this  CPT  location.  Figures  14  and  16  show  nearly  mirror  symmetrical  plots 
which  can  lead  to  non-unique  data  and  inversion  problems.  Non-unique  data  were 
encountered  in  some  of  the  simple,  early  developmental  models  which  used  only  a  few 
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symmetrical  rays  for  each  source  and  suggest  that  non-unique  data  can  arise  from  ray 
path  simulation  through  theoretical  models.  These  plots  suggest  that  some  constraint 
may  be  required  during  inversion  and  is  further  discussed  below. 


Numerical  Model  Phase  Shift  vs 
Spatially  Weighted  Straight  Ray  Phase  Shift 
30  Second  CPT  Period 

0-0.305  m  (0-1  ft)  Theoretical  Source  Location 


Phase  Shift 

0.6  0.5  0.4  0.3  0.2  0.1  0 


Numerical  Phase  Shift  from  FitAmpPhase  — ■ — Straight  Ray  Approximation  from  HydraulicTomAnal  •  Source  Location 


Figure  14  -  A  comparison  of  30-sec  CPT  period  phase  shift  values  from  a  numerical 
model  and  the  spatially  weighted  ray  path  method  at  the  0-0.305  m  (0-1  ft)  source 
location. 
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Numerical  Model  Phase  Shift  vs 
Spatially  Weighted  Straight  Ray  Phase  Shift 
30  Second  CPT  Period 

5.19-5.49  m  (17-18  ft)  Theoretical  Source  Location 


Phase  Shift 
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♦  Numerical  Phase  Shift  from  FitAmpPhase 


— ■ - Straight  Ray  Approximation  from  HydraulicTomAnal 


Source  Location 


Figure  15  -  A  comparison  of  30-sec  CPT  period  phase  shift  values  from  a  numerical 
model  and  the  spatially  weighted  ray  path  method  at  the  5.19-5.49  m  (17-18  ft)  source 
location. 


Numerical  Model  Phase  Shift  vs 
Spatially  Weighted  Straight  Ray  Phase  Shift 
30  Second  CPT  Period 

10.37-10.68  m  (34-35  ft)  Theoretical  Source  Location 


Phase  Shift 
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■  Straight  Ray  Approximation  from  HydraulicTomAnal  —♦^Numerical  Phase  Shift  from  FitAmpPhase  •  Source  Location 


Figure  16  -  A  comparison  of  30-sec  CPT  period  phase  shift  values  from  a  numerical 
model  and  the  spatially  weighted  ray  path  method  at  the  10.37-10.68  m  (34-35  ft)  source 
location. 
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After  the  goodness  of  fit  between  the  straight  ray  and  numerical  phase  shift  data 

were  evaluated  (Figs.  11-16),  the  straight  ray  phase  shift  data  (Table  3),  along  with  the 

numerical  phase  shift  data  (Table  4)  were  inverted  through  the  element  matrix  by  SVD 

analysis,  a  method  of  least  squares  fitting  and  inversion.  As  expected,  direct  inversion  of 

the  straight  ray  model  data  reproduced  the  input  model  K  values  for  each  of  the  layers 

with  practically  no  error  (Table  3).  The  percent  standard  deviation  on  the  K  values  for 

each  of  the  elements  were  essentially  zero,  implying  the  inversion  was  almost  perfect  for 

a  data  set  with  no  noise.  Random  error  of  5%  and  10%  was  applied  by  Monte  Carlo 

simulation  to  replicate  a  normal  and  worst-case  scenario  of  ambient  noise.  The  2.5%  - 

5.2%  range  indicates  the  inherent  error  associated  with  levels  of  random  noise  in  the 

middle  layer  and  a  4-sec  CPT  period  (Table  3).  In  contrast,  the  6.9%  -  14.5%  range  is  the 

inherent  error  associated  with  the  levels  of  random  noise  in  the  middle  layer  and  a  30-sec 

period  and  indicates  that  the  period  difference  tends  to  amplify  the  effect  of  random  error. 

Table  3  -  SVD  analysis  of  spatially  weighted  straight  ray  approximation  phase  shift 
through  a  three-element,  eight-node,  10.68  m  (35  ft)  thick  model  used  to  verify  the 
heterogeneous  extension. 


Spatially  Weighted  Straight  Ray  SVD  Analysis 

4-Sec  CPT  Period 

Monte  Carlo  -  No  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003 

3.616E-19 

0.00 

2 

0.006 

3.137E-18 

0.00 

3 

0.003 

3.870E-19 

0.00 

Monte  Carlo  -  5%  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003 

1.849E-05 

0.62 

2 

0.005999 

1.515E-04 

2.52 

3 

0.003 

1.899E-05 

0.63 

Monte  Carlo  -  1 0%  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.002999 

3.482E-05 

1.16 

2 

0.006018 

3.098E-04 

5.15 

3 

0.003001 

3.702E-05 

1.23 

Spatially  Weighted  Straight  Ray  SVD  Analysis 

30-Sec  CPT  Period 

Monte  Carlo  -  No  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003 

8.636E-19 

0.00 

2 

0.006 

7.492E-18 

0.00 

3 

0.003 

9.242E-19 

0.00 

Monte  Carlo  -  5%  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003 

5.063E-05 

1.69 

2 

0.006012 

4.159E-04 

6.92 

3 

0.003001 

5.205E-05 

1.73 

Monte  Carlo  - 10%  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003 

9.542E-05 

3.18 

2 

0.006109 

8.869E-04 

14.52 

3 

0.003005 

1.017E-04 

3.38 
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Inversion  of  the  4-sec  period  numerical  phase  shift  data  through  the  element 
matrix  was  also  reasonable  and  the  percent  standard  deviation  on  the  K  values  for  the 
middle  elements  was  3.1%  (Table  4),  in  the  absence  of  random  noise.  The  error 
associated  with  the  straight  ray  method  is  about  14.1  %  error  in  the  recovery  of  the  0.006 
ft/sec  K  by  the  straight  ray  method  (i.e.,  0.006  vs.  0.005s  ft/sec).  These  error  percents 
indicate  that  the  spatially  weighted  straight  ray  model  and  4-sec  CPT  period  can  resolve 
layers  of  about  1  m  (3  ft)  in  thickness  with  about  16-19%  total  error. 

Inversion  of  the  30-sec  period  numerical  phase  shift  data  through  the  element 
matrix  had  3.6%  percent  standard  deviation  on  the  K  values  for  the  middle  element 
(Table  4),  in  the  absence  of  random  noise.  This  inversion  was  constrained  slightly;  the 
offset  curves  (Fig.  15)  and  nearly  mirror  symmetric  plots  in  the  upper  and  lower  elements 
(Fig.  14  and  Fig.  16)  tend  to  suggest  non-uniqueness  data  issues  were  arising  during 
inversion.  The  SVD  analysis  was  slightly  weighted  with  a  constrained  least  squares 
factor  of  0.25,  which  gives  a  small  weight  to  the  initial  estimates  of  K  to  overcome  non¬ 
unique  data  and  shouldn’t  unnecessarily  restrain  the  analysis.  The  error  associated  with 
the  straight  ray  method  is  about  25%  error  in  the  recovery  of  the  0.006  ft/sec  K  by  the 
straight  ray  method  (i.e.,  0.006  vs.  0.0045  ft/sec).  These  error  percents  indicate  that  the 
spatially  weighted  straight  ray  model  and  30-sec  CPT  period  can  resolve  layers  of  about  1 
m  (3  ft)  in  thickness  with  about  27-29%  total  error. 
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Table  4  -  SVD  analysis  of  spatially  numerical  phase  shift  through  a  three-element,  eight- 
node,  10.68  m  (35  ft)  thick  model  used  to  verify  the  heterogeneous  extension. 


Numerical/Straight  Ray  Model  SVD  Analysis 

4  Sec  CPT  Period 

Monte  Carlo  No  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003032 

2.365E-05 

0.78 

2 

0.005155 

1.608E-04 

3.12 

3 

0.002986 

2.474E-05 

0.83 

Monte  Carlo  5%  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003032 

1.875E-05 

0.62 

2 

0.005154 

1.206E-04 

2.34 

3 

0.002986 

1.885E-05 

0.63 

Monte  Carlo  10%  Error 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.003031 

3.545E-05 

1.17 

2 

0.005169 

2.464E-04 

4.77 

3 

0.002987 

3.675E-05 

1.23 

Numerical/Straight  Ray  Model  SVD  Analysis 

30  Sec  CPT  Period 

Monte  Carlo  No  Error  -  CLS  0.25 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.002869 

4.125E-05 

1.44 

2 

0.004514 

1 .644E-04 

3.64 

3 

0.002847 

4.256E-05 

1.49 

Monte  Carlo  5%  Error  -  CLS  0.25 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.002869 

3.830E-05 

1.33 

2 

0.004513 

9.347E-05 

2.07 

3 

0.002848 

3.703E-05 

1.30 

Monte  Carlo  10%  Error  -  CLS  0.25 

Element 

K  (ft/sec) 

Stdr  Dev  K 

%  Stdr  Dev  K 

1 

0.00287 

7.678E-05 

2.67 

2 

0.004515 

1 .869E-04 

4.14 

3 

0.002849 

7.419E-05 

2.60 

Horizontal  Ray  Paths  -  ZOP  Profiles 

Introduction 

The  first  phase  of  this  project  was  to  use  horizontal  rays  only  and  thus  collect 
ZOP  profiles  where  the  source  and  receiver  were  at  the  same  elevations  in  the  source  and 
receiver  wells.  The  source  signal  was  generated  by  the  pneumatic  method,  as  described 
earlier,  using  a  manually  tuned  frequency  generator  giving  a  period  of  between  3-4 
seconds.  The  signal  frequency  was  manually  adjusted  to  give  a  signal  that  best 
represented  a  sinusoidal  form.  It  was  found  that  it  was  necessary  to  stay  near  the  natural 
frequency  of  the  well  for  best  results.  Two  kinds  of  profiles  for  the  receiver  well  were 
collected.  In  whole  well  tests  the  entire  column  of  the  source  well  was  oscillated  and 
only  the  receiver  well  location  was  packed  off.  Both  the  source  and  receiver  intervals 
were  isolated  by  straddle  packers  for  the  point  source  well  profdes.  Early  research  results 
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in  this  project  showed  that  it  was  necessary  to  pack  off  the  receiver  interval  to  obtain  the 
best  aquifer  response.  Both  the  stressed  interval  of  the  source  well  and  the  isolated 
receiver  interval  in  the  receiver  well  were  about  0.5  m  in  length.  The  locations  below  top 
of  casing  (BTOC)  were  referenced  to  the  center  of  the  stressed  or  received  interval.  Each 
location  center  was  approximately  0.3  m  from  the  next,  so  that  one  location  overlapped 
with  the  adjacent  locations.  The  overlapping  intervals  acted  much  like  a  centered  moving 
average,  where  the  vertical  changes  in  aquifer  heterogeneity  were  averaged  over  the  0.5 
m  interval,  but  were  assigned  to  the  center  point.  At  this  stage  we  were  developing 
processing  techniques  and  much  of  the  work  was  done  by  hand  in  multiple  steps.  The 
details  of  data  processing  are  given  in  Engard  et  al.  (2005)  and  Engard  (2006).  Later  this 
would  be  automated  in  the  FitAmpPhase  program,  which  would  speed  things  up 
considerably. 

One  is  able  to  approximate  the  diffusivity  from  the  final  corrected  amplitude 
derived  exponential  decay  term  d  and  the  phase  shift  Or,  equation  (18).  In  theory,  we 
have  two  independent  measures  of  K,  one  from  amplitude  and  one  from  phase 
measurements.  However,  we  have  found  the  amplitude  estimates  to  be  difficult  to  make 
because  we  do  not  know  precisely  the  effective  radius  to  use,  see  Engard  et  al.  (2005)  and 
Engard  (2006)  for  details.  For  this  reason,  we  will  only  present  here  results  of  K  for 
phase  measurements.  The  frequency  was  calculated  from  the  field  data  from  the 
reciprocal  of  the  fitted  source  well  period  for  each  CPT.  After  referring  to  the  literature 
an  initial  value  of  0.00001  was  used  for  Ss  (Fetter,  2001;  Domenico  and  Schwartz,  1998). 
Using  a  constant  value  for  Ss  is  unrealistic  but  is  necessary,  because  even  with  today’s 
technology,  it  is  difficult  to  measure  Ss  in  situ.  A  final  estimate  of  Ss  was  made  by 
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requiring  consistency  between  the  vertical  K  profiles  obtained  by  HRST  methods  and 
CPT  methods.  The  radial  distance  r  can  easily  be  measured  in  the  field  or  from  survey 
data.  With  some  algebraic  manipulation  estimates  of  K  can  be  made  from  the  CPT 
experimentally  measured  phase  and  amplitude  data.  Based  on  the  numerical  results 
presented  earlier,  the  CPT  derived  values  of  K  should  be  interpreted  as  distance  weighted 
averages  of  K  over  the  path  between  the  source  and  receiver  wells.  HRST  K  values  that 
differ  significantly  from  the  CPT  K  values  are  evidence  of  inter-well  heterogeneity. 

Results  From  High  Resolution  Slug  Testing  and  Continuous  Pulse  Testing 

For  this  project,  high-resolution  slug  test  (HRST)  techniques  (discussed  earlier) 
were  applied  to  newly  installed  wells  HT-1,  HT-2,  and  HT-3  after  they  were  properly 
developed.  HRST  data  from  other  wells  (Ross,  2004)  also  was  used  for  comparison  to 
continuous  pulse  tests  CPTs.  A  dual  packer  arrangement  with  a  0.5  m  interval  open  to 
the  formation  was  used  for  pneumatic  slug  testing  (Ross,  2004;  Zemansky  and  McElwee, 
2005;  Ross  and  McElwee,  2007). 

The  research  presented  here  uses  continuous,  controlled,  sinusoidal  pressure 
signals  as  a  means  to  estimate  vertical  profiles  of  well-to-well  hydraulic  diffusivity.  The 
received  signal  is  measured  at  various  depths  in  observation  wells  at  various  distances 
and  locations.  The  length  of  the  vertical  profiles  measured  by  the  CPT  methods  are 
limited  by  the  amount  of  open  screen  common  to  the  well  pair  in  question  and  by  the 
length  of  the  bottom  packer  on  the  source  and  receiver  double  packer  apparatus. 
Typically,  the  CPT  profile  was  about  8  m  in  length. 
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In  total,  7  line  source  (whole  well  tests,  the  entire  column  was  oscillated)  well 
pairs  were  tested  with  the  pneumatic  CPT  method  at  GEMS.  The  shortest  well  separation 
distance,  4.36  m,  was  between  well  HT-2  and  well  HT-3.  The  longest  separation 
distance,  1 1.5  m,  recorded  was  between  well  00-3  and  well  7-1.  These  results  are 
presented  in  figures  17-23.  The  averages  of  the  HRST  values  at  each  depth  for  the 
source  and  receiver  wells  are  plotted  along  with  the  highest  and  lowest  values  shown  by 
error  bars.  This  curve  is  labeled  HRST.  The  other  curve  labeled  CPT  presents  the  results 
of  analyzing  the  CPT  data  for  K.  The  two  curves  generally  agree  fairly  well,  with  the 
exception  of  figure  18.  It  appears  that  the  general  features  of  the  HRST  curve  are 
captured  by  the  CPT  curve,  but  it  seems  smoother.  This  is  probably  because  of  the  long 
line  source  geometry  giving  poorer  resolution.  It  is  unknown  at  this  time  why  the  CPT 
curve  in  figure  1 8  is  so  flat;  perhaps  it  is  due  to  some  experimental  or  processing  problem 
we  have  not  discovered.  In  general  it  is  difficult  to  analyze  the  whole  well  tests,  since  the 
source  area  is  so  large  and  multiple  paths  of  energy  lead  from  the  source  to  the  point 
receiver.  In  later  work  we  decided  not  to  use  whole  well  tests  due  to  the  difficulty  in 
analyzing  them. 
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Well  00-3  to  7-1 


0  0.001  0.002  0.003  0.004 

K  m/s 

Figure  17.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  line  source  pneumatic  CPT  between  wells  00-3  to  7-1. 
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Well  11-1  to  7-1 
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Figure  18.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profile  from  a  line  source  pneumatic  CPT  between  wells  00-3  to  7-1. 
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Well  HT-1  to  7-1 
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Figure  19.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  line  source  pneumatic  CPT  between  wells  HT-1  to  7-1. 
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Well  HT-2  to  7-1 
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Figure  20.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profile  from  a  line  source  pneumatic  CPT  between  wells  HT-2  to  7-1. 
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Well  7-1  to  HT-3 
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Figure  21 .  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  line  source  pneumatic  CPT  between  wells  7-1  to  HT-3. 
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Well  HT-1  to  HT-3 


—  CPT  -  HRST 
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Figure  22.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  line  source  pneumatic  CPT  between  wells  HT-1  to  HT-3. 
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Well  HT-2  to  HT-3 

-CPT-HRST 
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Figure  23.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  line  source  pneumatic  CPT  between  wells  HT-2  to  HT-3. 
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Five  point  source  profiles  (both  source  and  receiver  packed  off)  were  completed 
at  GEMS  with  the  pneumatic  CPT  method.  Also,  one  point  source  profile  was  completed 
with  an  early  crude  version  of  the  mechanical  pumping  or  injection  CPT  method  (an 
improved  version  of  the  mechanical  pumping  system  will  be  the  subject  of  later  sections 
of  this  report).  The  shortest  well  separation  distance  of  4.36  m  was  between  well  HT-2 
and  well  HT-3.  The  longest  separation  distance,  6.91  m,  was  recorded  between  well  HT- 
2  and  well  7-1.  These  results  are  presented  in  figures  24-  28  for  the  pneumatic  profiles 
and  in  figure  29  for  the  injection  profile.  The  presentation  style  is  the  same  as  for  figures 
17-23  with  the  HRST  curves  being  the  same  as  before.  The  CPT  curves  are  now  for  the 
point  source  CPT  method  and  seem  to  have  more  detail  and  are  more  closely  correlated 
to  the  HRST  data.  It  appears  that  the  point  source  CPT  tests  are  giving  better  K 
resolution,  as  we  might  expect.  Comparison  of  the  pneumatic  method  of  figure  24  and 
the  injection  method  of  figure  29  for  the  same  well  pair  shows  that  the  results  are  similar, 
but  some  differences  do  occur. 
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Well  11-1  to  7-1 
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Figure  24.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  point  source  pneumatic  CPT  between  wells  11-1  to  7-1. 


64 


Well  HT-1  to  HT-3 
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Figure  25.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profile  from  a  point  source  pneumatic  CPT  between  wells  HT-1  to  HT-3. 
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Well  HT-2  to  HT-3 
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Figure  26.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profde  from  a  point  source  pneumatic  CPT  between  wells  HT-2  to  HT-3. 
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Well  HT-2  to  7-1 
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Figure  27.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profile  from  a  point  source  pneumatic  CPT  between  wells  HT-2  to  7-1. 
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Well  7-1  to  HT-3 


-CPT-HRST 


0  0.001  0.002  0.003  0.004 

K  m/s 


Figure  28.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profile  from  a  point  source  pneumatic  CPT  between  wells  7-1  to  HT-3. 
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Figure  29.  Comparisons  of  the  composite  HRST  values  from  each  well  and  the  average 
estimated  K  profile  from  a  point  source  injection  CPT  between  wells  11-1  to  7-13. 


It  is  evident  that  the  CPT  profiles  mimic  the  general  trends  in  the  HRST  K 
profiles  measured  at  the  respective  wells.  Overall,  the  CPT  data  appear  to  average  the  K 
profiles  of  the  well  pair  in  question.  However,  there  are  important  differences.  The 
heterogeneities  of  the  geologic  material  between  the  well  pair  are  probably  the  cause  of 
this  difference;  and  the  difference  can  not  be  fully  explained  without  using  more 
advanced  models  and  numerical  solutions.  The  point  source  data  appear  to  increase  the 
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resolution  of  the  data,  distinguishing  variations  in  K  that  are  not  present  in  the  HRST  and 
line  source  data. 

Calculation  of  “Anomalous  K  Values” 

The  previously  discussed  comparison  of  the  straight  ray  approximation  to 
numerical  results  indicates  that  using  a  spatially  weighted  average  for  K  should  be 
appropriate.  This  allows  the  interpretation  of  the  hydraulic  conductivity  (K)  data  that 
have  been  collected  for  high-resolution  slug  tests  at  wells  and  the  inter- well  spatially 
averaged  K  that  has  been  determined  by  the  continuous  pulse  testing  horizontal  ray  path 
data.  Unfortunately,  there  is  no  unique  way  to  do  the  spatial  weighting  with  the 
horizontal  ray  path  data,  since  we  only  have  one  ray  path  crossing  each  segment  of  the 
aquifer.  When  we  are  able  to  collect  diagonal  ray  path  data,  we  will  have  multiple  rays 
that  cross  each  given  segment  of  the  aquifer  and  it  will  be  possible  to  estimate  the  spatial 
averaging  to  some  scale  limited  by  the  density  of  ray  paths.  However,  in  the  present  case 
(horizontal  ray  path  data)  one  must  assume  some  spatial  averaging  scheme  to  interpret 
what  the  inter-well  average  K  is  telling  us  about  the  variation  of  K  between  wells.  It  is 
well  known  that  slug  tests  only  give  a  K  value  that  is  representative  near  the  well. 
Therefore,  we  should  give  less  weight  to  the  slug  test  values  and  more  weight  to  the  inter¬ 
well  average  K  determined  from  the  horizontal  ray-path  data.  Arbitrarily,  let  us  assume 
that  the  weight  for  each  of  the  two  slug  test  values  is  1/6  and  the  weight  for  the  inter- well 
average  K  is  2/3  =  4/6.  These  weighting  coefficients  add  up  to  1,  as  any  weighting 
scheme  should.  This  assumption  will  allow  us  to  calculate  a  new  value  of  K  between  the 
source  and  receiver  wells  that  may  be  different  from  the  slug  test  K  values  or  the  inter- 
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well  average  K  determined  from  the  horizontal  ray  path  data.  In  what  follows  we  will 
call  this  calculated  value  of  K  between  the  source  and  receiver  wells  the  “anomalous  K 
value”  (Engard  et  al.,  2006)  because  it  can  be  different  from  any  of  the  experimentally 
determined  K  values.  Ideally,  if  the  K  values  changed  in  a  linear  fashion  from  the  source 
well  to  the  receiver  well,  the  inter-well  average  K  determined  from  the  horizontal  ray 
path  data  should  fall  between  the  values  of  K  determined  by  slug  tests  at  each  well. 
However,  we  have  observed  that  this  often  is  not  the  case,  which  means  that  the  K  values 
are  not  varying  linearly  between  wells.  The  above  outlined  scheme  allows  us  to  calculate 
an  “anomalous  K  value”  that  shows  how  K  may  in  fact  be  varying  between  wells. 

This  procedure  has  been  applied  to  three  of  the  source-receiver  well  pairs  from 
which  we  have  collected  data.  The  calculations  for  the  Well  pairs  HT-1  to  HT-3,  HT-2 
to  HT-3,  and  7-1  to  HT-3  are  show  in  figures  30,  31,  and  32  respectively.  What  we 
observe  is  that  the  inter-well  average  K  (CPT  K  value  in  the  figures,  shown  in  dark  blue) 
is  many  times  outside  the  interval  defined  by  the  two  slug  test  values  of  K  shown  in  pink 
and  yellow.  This  means  that  the  K  value  is  not  varying  linearly  between  wells  and  an 
anomalous  value  outside  that  range  may  be  calculated.  The  anomalous  values  are  shown 
in  light  blue-green.  When  the  anomalous  curve  is  significantly  outside  the  slug  test  K 
interval,  we  have  an  indication  that  the  K  value  between  wells  is  significantly  different 
from  that  observed  at  the  source  and  receiver  wells.  The  values  of  the  anomalous  K  are 
calculated  using  the  weighting  scheme  detailed  earlier.  It  should  be  reiterated  that  the 
choice  of  weighting  scheme  is  arbitrary  at  this  point  and  is  not  unique.  However,  the 
calculations  presented  here  should  be  useful  in  identifying  areas  of  “anomalous  K” 
between  the  source  and  receiver  wells. 
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Calculated  Anomalous  K  Values 
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Figure  30.  Experimentally  determined  K  values  from  high  resolution  slug  tests  and 
horizontal  ray  path  data  along  with  “anomalous  K  values”  calculated  from  the  weighting 
scheme,  for  well  pair  HT-1  and  HT-3. 


Calculated  Anomalous  K  Values 
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Figure  3 1 .  Experimentally  determined  K  values  from  high  resolution  slug  tests  and 
horizontal  ray  path  data  along  with  “anomalous  K  values”  calculated  from  the  weighting 
scheme,  for  well  pair  HT-2  and  HT-3. 
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Calculated  Anomalous  K  Values 
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Figure  32.  Experimentally  determined  K  values  from  high  resolution  slug  tests  and 
horizontal  ray  path  data  along  with  “anomalous  K  values”  calculated  from  the  weighting 
scheme,  for  well  pair  7-1  and  HT-3. 

Reproducibility  and  Reciprocity 

We  have  investigated  the  reproducibility  of  the  data  and  the  reciprocity  of  source 
and  receiver  wells  (Engard  et  al.,  2006).  In  well  pair  HT-1  and  HT-3  we  have  taken  data 
at  two  different  times  and  with  the  source  and  receiver  locations  reversed.  The  results  are 
shown  in  Figure  33  for  the  measured  phase  which  is  the  basic  data  used  to  calculate  K.  It 
is  seen  that  the  signals  are  reproducible  within  experimental  error  over  an  extended  time 
interval  between  data  collection  and  with  the  source  and  receiver  reversed.  Similarly,  for 
well  pair  HT-2  and  HT-3  we  have  taken  data  at  four  different  times  and  once  with  the 
source  and  receiver  reversed.  These  data  sets  are  shown  in  Figure  34.  The  general  shape 
and  features  of  the  phase  shift  curve  are  reproduced  and  the  various  data  sets  agree  within 
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experimental  accuracy.  Therefore,  we  conclude  that  the  data  are  reproducible  and  that 
they  are  nearly  independent  of  source  and  receiver  position,  within  experimental  error. 


Figure  33.  Phase  data  for  well  pair  HT-1  and  HT-3  at  two  different  times  with  the  source 
and  receiver  reversed. 
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Figure  34.  Phase  data  for  well  pair  HT-2  and  HT-3  at  four  different  times  with  the  source 
and  receiver  reversed  on  one  data  set. 


Diagonal  Ray  Paths  -3-4  Sec.  Pneumatic  MOG  Data 

Introduction 

In  the  second  phase  of  this  research  we  started  collecting  MOG  data  with  diagonal 
ray  path  using  the  pneumatic  method  described  earlier.  It  is  only  with  MOG  data  that  a 
true  tomography  can  be  utilized,  since  it  is  necessary  to  have  multiple  rays  at  different 
angles  passing  through  the  area  of  interest.  The  period  of  the  oscillations  was  determined 
by  a  manually  set  frequency  generator,  so  the  period  was  not  absolutely  constant  between 
records,  but  was  in  the  3-4  sec.  range.  At  processing  time  all  the  records  were  corrected 
to  a  standard  period.  Figure  8  shows  all  the  hydraulic  tomography  (HT)  wells  installed  as 
part  of  this  project.  HT-3  sits  at  the  center  of  the  array  with  other  wells  surrounding  it. 
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All  pairs  of  the  hydraulic  tomography  wells  using  HT-3  as  a  common  center  were  tested 
using  the  pneumatic  method.  This  allows  a  series  of  slices  through  the  aquifer,  which 
when  taken  as  a  whole  allows  a  3-D  view  of  the  hydraulic  conductivity  distribution  as 
determined  by  this  research.  The  following  sections  summarize  the  model  selection  for 
inversion  and  the  results  for  constrained  inversion  of  the  data.  For  a  more  complete 
description  of  the  details  see  McElwee  et  al.,  2007;  Wachter,  2008;  Wachter  et  al.,  2008. 

Sensitivity  of  Various  Models  to  Inversion  for  K  Values 

It  is  well  known  that  the  stability  and  uniqueness  of  an  inverse  problem  depends 
both  on  the  data  collected  (quantity  and  quality)  and  the  structure  of  the  selected  model 
for  inversion.  The  number  of  parameters  to  be  estimated  and  the  size  and  shape  of  the 
zone  they  occupy  are  critical  to  the  inversion.  For  these  reasons  a  number  of  runs  were 
made  with  synthetic  and  real  data  to  try  and  determine  a  reasonable  model  structure  for 
our  inversion  (Wachter,  2008). 

Theoretical  values  of  phase  and  amplitude  for  more  complex  models  were  run 
through  data  processing  programs  before  applying  the  programs  to  field  data.  The 
synthetic  data  sets  could  be  generated  with  no  error  or  a  given  level  of  random  noise. 

The  current  version  of  the  SVD  inverse  program  has  the  ability  to  perform  Monte  Carlo 
simulations  in  the  presence  of  random  error.  The  Monte  Carlo  simulations  were  run  with 
both  +/-  1%  and  +/-  5%  noise  for  1000  simulations.  A  variety  of  models  with  differing 
numbers  of  nodes  in  the  x  and  z  directions,  numbers  of  elements,  and  numbers  of  zones 
were  investigated.  All  initial  models  used  100  ray  paths,  consisting  of  10  source 
locations  and  10  receiver  locations,  and  Ss  values  of  0.00018  or  .00001.  K  values  were 
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arbitrarily  chosen  to  start  at  0.000914  m/s  (0.003  ft/s)  at  shallow  depths  and  gradually 
increased  with  depth  to  0.00213  m/s  (0.007  ft/s).  Although  the  K  values  were  arbitrarily 
chosen  for  the  modeling  phase,  they  fall  within  the  range  observed  at  the  site  (0.000305 
m/s  to  0.00305  m/s)  from  HRST  and  other  methods. 

The  HydTomAnal  program  calculates  the  length  of  each  ray  path  through  a 
particular  element  or  near  a  given  node.  The  total  amount  of  ray  path  length  through 
each  element  or  associated  with  each  node  can  be  calculated  by  adding  the  lengths  from 
each  individual  ray  path  of  the  100  rays  used.  The  ray  path  sums  give  a  measure  of  the 
sensitivity  of  a  given  model  to  the  K  value  in  an  element  or  near  a  node.  The  ray  path 
density  was  highest  in  the  center  of  the  region,  so  there  was  less  resolution  at  the  top  and 
bottom  of  the  modeled  area.  The  problem  can  be  avoided  by  having  spatially  variable 
element  sizes  across  the  model.  The  latest  versions  of  the  processing  programs  offer  the 
ability  to  specify  K  by  zones,  which  are  formed  by  one  or  more  nodes  or  elements  and 
must  be  input  manually.  The  purpose  of  the  zones  is  to  provide  variable  resolution  across 
the  model,  with  finer  zones  towards  the  center  and  coarser  zones  at  the  edges  of  the  grid 
where  fewer  rays  are  crossing. 

The  amount  of  error  produced  by  a  given  set  of  input  parameters  was  balanced 
with  the  amount  of  resolution  provided  by  that  particular  model.  Increasing  the  number 
of  zones  increases  the  resolution,  but  only  at  the  cost  of  increased  error.  A  mode  with  50 
elements  and  1 6  zones  resulted  in  the  least  amount  of  error  of  the  models  studied  here, 
but  some  models  with  more  zones  also  produced  acceptable  amounts  of  error.  A  good 
balance  of  error  and  resolution  when  Ss  equaled  0.00001  was  achieved  with  this  16  zone, 
50  element  model  (6  nodes  in  x  direction).  The  average  error  was  7.79%  in  the  presence 
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of  1%  noise.  This  16  zone  50  element  model  with  six  nodes  in  the  x  direction  will  be 


used  in  all  following  discussions. 

The  grid  layout  for  the  chosen  16  zone  50  element  model  is  shown  in  Figure  35. 
In  Figure  35,  each  of  the  50  elements  is  numbered,  with  element  1  at  the  bottom  of  the 
source  well  and  element  50  at  the  top  of  the  receiver  well.  The  greatest  resolution  is 
provided  in  the  middle  of  the  grid  while  the  top  and  bottom  have  the  least  resolution.  As 
discussed  above  the  total  ray  path  length  associated  with  each  element  or  zone  is  a 
measure  of  the  model  sensitivity  to  the  value  of  K  in  that  element  or  zone.  The  sum  of 
ray  path  lengths  going  through  each  zone  of  the  chosen  model  for  the  suite  of  rays  used 
was  also  calculated  and  is  presented  below  in  Figure  36.  They  were  calculated  using 
field  geometry  and  the  actual  number  of  rays  collected  in  the  field  for  each  well  pair. 
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Figure  35:  The  grid  shows  the  division  of  elements  into  16  zones. 


In  each  well  pair,  the  center  of  the  model  prior  to  zoning  had  the  highest  sum  of 
ray  path  lengths  because  the  most  ray  paths  passed  through  those  areas.  Other  elements 
were  combined  together  to  produce  zones  with  sums  comparable  to  the  value  in  the 
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center.  Zone  sums  may  differ  by  a  factor  of  two  but  should  not  vary  by  as  much  as  an 
order  of  magnitude  for  a  well  behaved  model.  Zone  sums  in  a  given  zone  are  fairly 
similar  between  well  pairs  with  approximately  the  same  number  of  ray  paths.  Variation 
occurs  from  one  well  pair  to  another  partially  due  to  differing  radii  but  largely  due  to  the 
changing  number  of  ray  paths. 
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(d).  HT-5  to  HT-3  (190  rays) 
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(e).  HT-6  to  HT-3  (300  rays) 
361.10 

454.76 
143.64 
202.04 
210.42 
158.90 


214.51 

212,24 

209.94 

211.12 


209.34 

211.15 

211.16 
216.29 


Figure  36:  The  sums  of  ray  paths  in  each  zone  for  well  pairs  HT-3  to  HT-2  (a),  HT-3  to 
HT-3  (b),  HT-4  to  HT-3  (c),  HT-5  to  HT-3  (d),  and  HT-6  to  HT-3  (e). 


SVD  Processing 

After  choosing  a  zoning  model,  field  data  were  run  through  the  inversion  program 
to  determine  K  values.  As  stated  in  the  Modeling  section,  the  model  chosen  for  this 
scenario  was  a  16  zone  model  where  K  was  determined  by  elements  and  there  were  six 
nodes  in  the  X  direction.  The  data  set  from  well  HT-3  to  well  HT-2  was  first  examined 
because  it  seemed  to  be  the  best  of  the  initial  data  sets.  The  HRST  K  values  determined 
in  previous  tests  were  input  as  constant  K  nodes  to  help  fix  the  other  K  values  within  a 
reasonable  range.  HRST  results  were  processed  by  Brett  Engard  (2006)  for  wells  HT-1, 
HT-2,  and  HT-3,  and  by  Pema  Deki  (2008)  for  wells  HT-4,  HT-5,  and  HT-6  (Appendix 
B). 

Contour  plots  were  made  of  K  values  plotted  against  elevation  and  the  radial 
distance  between  wells  using  a  program  called  QuikGrid,  a  public  domain  program.  The 
program  contours  between  points  written  in  an  x,y,z  format,  in  this  case  corresponding  to 
radius,  elevation,  and  K.  The  contour  interval  chosen  is  0.0002  m/s.  The  HRST  values 
were  chosen  for  the  K  values  at  each  well  in  the  plot.  Interwell  K  values  were 
determined  by  the  SVD  analysis,  a  method  using  least  squares.  In  the  contour  plots  of  K, 
the  source  well  is  on  the  left  side  and  the  receiver  well  is  on  the  right  side. 

Initially  the  field  data  were  processed  using  an  unconstrained  SVD  procedure. 

The  results  were  unstable  with  regions  of  K  occurring  that  were  known  to  be 
unreasonable.  The  SVD  inverse  program  performs  perfectly  on  model  data  without 
noise,  so  it  must  be  much  more  sensitive  to  noise  than  originally  thought.  To  compensate 
for  the  sensitivity  to  noise,  a  seven  point  filter  was  used  on  the  data.  In  addition,  noise 
reduction  was  attempted  by  editing  larger  offset  rays,  where  noise  was  expected  to  be 
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greater.  Little  improvement  was  observed  due  to  filtering  and  ray  path  editing. 
Apparently,  the  inverse  procedure  needed  some  additional  conditioning  to  become  stable. 
As  an  alternative  processing  scheme,  an  SVD  least  squares  fit  was  employed  constrained 
by  the  HRST  data,  which  is  detailed  in  the  following  section. 

Constrained  SVD  Results 

Inverse  problems  are  commonly  constrained  with  data  known  from  other  sources 
or  methods;  in  this  case,  HRST  results  were  used  to  constrain  the  inversion  for  K  values. 
Initial  guesses  of  K  in  each  zone  were  obtained  through  a  linear  interpolation  of  HRST 
values  at  the  same  Z  elevations.  The  sum  of  squared  errors  (SSE)  was  calculated  by 
comparing  the  phase  values  measured  in  the  field  to  the  phase  values  calculated  using 
SVD.  A  weighting  factor  is  used  in  the  latest  version  of  the  SVD  program  to  determine 
to  what  extent  the  HRST  results  constrain  the  inversion.  A  weighting  factor  of  zero  is 
equivalent  to  the  unconstrained  SVD  analysis,  and  increasing  values  for  the  weighting 
factor  result  in  increasing  weight  given  to  the  HRST  results  and  therefore  less  deviation 
from  HRST  values.  For  this  study,  a  factor  of  1.0  was  used,  resulting  in  about  equal 
weight  of  the  HRST  data  and  the  tomographic  data.  The  K  value  in  a  zone  is  only 
changed  if  it  is  still  in  the  approximate  range  of  values  seen  from  HRST.  The  results 
were  calculated  using  two  values  for  Ss,  10"5  (Table  5)  and  1.5  x  10"5  (Table  6). 
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Contour  plots  were  made  of  K  values  plotted  against  elevation  and  the  radial 
distance  between  wells  for  the  data  constrained  by  the  HRST  results.  The  HRST  values 
are  used  at  the  left  and  right  ends  of  the  plot,  with  the  source  on  the  left  and  the  receiver 
on  the  right.  Interwell  K  values  in  the  following  plots  were  all  determined  by  the 
constrained  SVD  analysis.  The  phase  is  a  ratio  between  Ss  and  K,  so  changes  in  Ss  will 
also  result  in  changes  in  K.  This  introduces  a  potential  source  of  error  because,  due  to  the 
difficulty  of  measuring  Ss  in  situ,  a  value  was  obtained  from  the  literature  rather  than 
from  field  measurements.  To  investigate  the  effect  of  Ss,  the  constrained  SVD  analysis 
was  conducted  on  all  of  the  data  using  Ss  values  of  both  10"5  (Figures  37-44)  and  1.5  x 
10~5  (Figures  45-52).  A  value  of  1.5  x  10~5  in  general  results  in  smoother  transitions 
between  zones.  The  negative  aspect  of  choosing  the  higher  Ss  value  is  that  well  pair  HT- 
6  to  HT-3,  which  already  had  higher  than  expected  K  values  with  the  lower  Ss  (Figure 
44),  continues  to  increase  above  the  expected  range  (Figure  52). 

Based  on  other  work  at  the  site,  and  in  particular  HRST,  K  values  at  the  site  are 
known  to  range  from  approximately  0.0003  m/s  up  to  0.003  m/s.  The  K  values  in  figures 
37  and  45  are  all  within  this  range.  The  trend  also  matches  that  seen  in  HRST  results, 
with  low  K  material  near  the  top,  a  high  K  region  in  the  middle,  another  high  K  region 
beginning  at  the  bottom  of  the  plot,  and  a  low  K  zone  between  the  two  high  K  regions. 
The  data  set  from  HT-3  to  HT-2  was  used  to  verify  that  the  program  was  working 
correctly  before  extending  the  analysis  to  other  well  pairs.  Figure  45,  using  a  value  of  1.5 
x  10'5  for  Ss,  shows  a  smoother  transition  between  points  than  Figure  37,  which  is 
physically  a  more  likely  scenario. 
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Plots  were  also  made  of  the  HT-3  to  HT-2  data  set  using  less  than  750  rays  to 
determine  if  fewer  rays  can  provide  the  same  results.  The  number  of  rays  in  each 
example  was  based  on  the  ray  path  geometry  of  the  other  well  pairs.  The  270  ray  path 
example  used  all  receiver  data  for  each  source  used  but  only  every  third  source  location, 
just  like  well  pair  HT-6  to  HT-3.  Similarly,  the  170  ray  path  example  followed  the  same 
pattern  as  well  pair  HT-5  to  HT-3  and  the  90  ray  path  example  followed  the  pattern  of 
well  pair  HT-4  to  HT-3.  The  three  following  figures  (Figures  38-40)  show  the  same 
trend  seen  in  Figure  37,  but  the  magnitudes  of  the  K  values  decrease  as  the  number  of  ray 
paths  decreases.  The  270  ray  path  scenario  (Figure  38)  is  closest  to  the  750  ray  path 
scenario.  The  bottom  zone  is  about  the  same  in  the  750  and  270  ray  path  cases,  but  the  K 
values  in  the  bottom  zone  are  noticeably  smaller  in  the  two  cases  with  less  ray  paths.  The 
plots  using  the  higher  Ss  value  (Figures  46-48)  also  show  the  same  trends,  but  the 
transitions  between  zones  are  smoother. 

The  data  set  presented  in  Figure  41  is  not  as  accurate  as  the  other  data  sets.  The 
amount  of  error  between  calculated  and  observed  phases  was  greater  than  that  in  other 
well  pairs.  Problems  with  this  data  set  are  likely  caused  by  the  nitrogen  leaks  at  the  time 
of  data  collection.  The  equipment  was  repaired  after  this  data  set  was  completed.  In 
spite  of  the  problems,  the  plot  shows  the  same  general  zones  of  high  and  low  K  seen 
elsewhere.  Figure  49,  using  the  larger  Ss  value,  depicts  the  same  zones  of  high  and  low 
K.  The  processing  program  is  probably  not  causing  the  problems  because  it  has  been 
constrained  and  other  well  pairs  do  not  have  as  many  problems.  So,  drawing  definite 
conclusions  about  this  well  pair  would  likely  require  recollecting  the  data  with  the  current 
repaired  equipment. 
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Both  of  the  original  well  pairs  verified  the  use  of  the  constrained  processing 
program  by  showing  the  trends  observed  in  HRST  results,  so  the  data  from  fall  2007  were 
examined  for  the  three  newest  wells.  The  vertical  intervals  were  varied  in  the  source  and 
receiver  wells  for  some  of  the  wells  pairs;  this  offered  the  opportunity  to  determine  if 
data  were  being  collected  at  adequate  spatial  intervals  for  appropriate  resolution.  The 
data  in  figures  42  and  50  show  the  overall  trend  observed  between  HT-3  and  HT-2  and 
between  HT-3  and  HT-1.  Once  again,  the  plots  demonstrate  the  expected  trends  of  high 
and  low  K  zones.  The  high  K  zone  near  the  top  of  the  plot  is  not  seen  elsewhere  in  that 
portion  of  the  aquifer,  but  the  values  are  at  least  within  the  overall  range  determined  by 
other  methods.  This  could  potentially  be  caused  by  a  combination  of  previously 
discussed  problems  of  resolution  in  the  top  of  the  sampling  area  in  combination  with  the 
low  number  of  ray  paths  used  for  this  particular  well  pair  (100,  compared  to  750  for  HT-3 
to  HT-2). 

The  results  between  well  HT-5  and  well  HT-3  are  presented  in  Figures  43  and  51. 
The  difference  between  the  two  figures  is  that  the  transitions  between  K  values  are 
smoother  in  the  plot  using  the  higher  value  for  Ss.  Some  of  the  values  at  the  bottom  of 
the  plot  are  slightly  above  the  general  expected  range,  but  still  within  reason.  K  values 
have  been  shown  to  slightly  exceed  0.003  m/s  in  some  of  the  HRST  data  toward  the 
bottom  of  the  wells.  The  same  trend  of  low  K  material  at  the  top,  a  moderately  high  K 
zone  in  the  middle,  and  high  K  material  at  the  bottom  is  again  observed  in  this  well  pair. 
As  with  the  plot  from  well  HT-4  to  well  HT-3,  the  relatively  large  region  of  very  low  K 
values  at  the  top  could  be  due  to  the  lower  number  of  ray  paths  for  this  well  pair. 
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The  contour  plots  from  well  HT-6  to  HT-3  (Figures  44  and  52)  show  the  same 
trend  seen  in  tomography  experiments  between  all  the  other  well  pairs,  as  well  as  in 
HRST  results.  The  zone  of  very  high  K  in  the  middle  left  side  of  the  plot  exceeds  the 
range  of  expected  values  for  the  site.  The  location  of  the  zone  could  be  due  to  the  survey 
design  for  this  well  pair,  namely,  the  source  locations  were  sampled  at  a  coarser  interval 
than  that  used  for  the  receiver  locations.  This  well  pair  is  the  only  one  examined  in  this 
study  for  which  increasing  Ss  resulted  in  K  values  farther  above  the  expected  range. 
Despite  the  problem  of  larger  than  expected  K  values,  the  transitions  between  zones  are 
again  smoother  using  1.5  x  10'5  instead  of  10'5  for  Ss. 

The  number  of  ray  paths  collected  for  a  well  pair  correlated  well  with  the 
reasonableness  of  the  K  values.  The  well  pair  with  the  best  results,  HT-3  to  HT-2  (750 
rays),  was  characterized  by  the  most  rays  of  any  of  the  well  pairs,  with  the  exception  of 
the  well  pair  with  equipment  problems.  The  well  pairs  of  HT-4  to  HT-3  and  HT-5  to  HT- 
3  had  100  rays  and  190  rays,  respectively,  and  some  of  the  higher  elevation  zones  were 
somewhat  lower  than  expected  for  the  site.  The  results  suggest  that  190  ray  paths  are  not 
enough  for  accurate  results.  Time  constraints  may  not  always  allow  for  750  ray  paths, 
but  there  does  seem  to  be  a  strong  correlation  with  the  accuracy  of  the  processing  results 
and  the  number  of  ray  paths.  The  work  with  editing  ray  paths  for  the  data  set  from  HT-3 
to  HT-2  also  lends  support  for  collecting  as  many  ray  paths  as  time  permits.  The 
resolution  of  K  values  decreased  as  more  rays  were  edited  out.  Although  it  takes  less 
time  to  collect  300  ray  paths  than  to  collect  750  ray  paths,  the  additional  rays  will  provide 
some  increase  in  accuracy. 
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Figure  37:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (750  rays,  Ss=  10'5). 
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Figure  38:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (270  rays,  Ss=  10~5). 
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Figure  39:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (170  rays,  Ss=  10'5). 
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Figure  40:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (90  rays,  Ss=  10'5). 
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Figure  41:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-1  as  the  receiver  well  (Ss=  10"5). 
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Figure  42:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-4  as  the  source 
well  and  HT-3  as  the  receiver  well  (Ss=  10"5). 
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Figure  43:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-5  as  the  source 
well  and  HT-3  as  the  receiver  well  (Ss=  10"5). 
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Figure  44:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-6  as  the  source 
well  and  HT-3  as  the  receiver  well  (Ss=  10"5). 
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Figure  45:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (750  rays,  Ss=  1.5  x  10'5). 
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Figure  46:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (270  rays,  Ss=  1.5  x  10'5). 
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Figure  47:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (170  rays,  Ss=  1.5  x  10'5). 
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Figure  48:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-2  as  the  receiver  well  (90  rays,  Ss=  1.5  x  10'5). 
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Figure  49:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-3  as  the  source 
well  and  HT-1  as  the  receiver  well  (Ss=  1.5  x  10'5). 
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Figure  50:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-4  as  the  source 
well  and  HT-3  as  the  receiver  well  (Ss=  1.5  x  10~5). 
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Figure  51:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-5  as  the  source 
well  and  HT-3  as  the  receiver  well  (Ss=  1.5  x  10~5). 
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Figure  52:  Interwell  K  values  from  constrained  SVD  analysis  with  HT-6  as  the  source 
well  and  HT-3  as  the  receiver  well  (Ss=  1.5  x  10"5). 
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Diagonal  Ray  Paths  -30  Sec.  Mechanical  Pumping  and  3  Sec.  Pneumatic 
Geoprobe  Source  MOG  Data 

Introduction 

The  last  phase  of  this  project  was  to  collect  data  with  the  mechanical  pumping 
system  and  to  see  if  the  sinusoidal  source  could  be  operated  from  the  end  of  a  direct  push 
system  such  as  Geoprobe.  Data  were  collected  using  the  same  set  of  hydraulic 
tomography  (HT)  wells  as  used  for  the  pneumatic  data  collection.  For  the  location  of  the 
Geoprobe  source  well  we  picked  a  location  that  keeps  HT-3  as  the  center  and  fills  in  a 
vacant  area,  as  shown  in  Figure  8.  The  pumping  system  was  constrained  to  a  period  of 
30  seconds  because  that  was  the  minimum  limit  of  our  computer  controlled  servo-valve. 

It  would  be  desirable  to  design  and  build  a  servo-valve  that  could  bridge  the  gap  between 
the  30  second  data  generated  here  and  the  3  second  data  of  the  pneumatic  system.  As 
discussed  earlier  the  wavelength  of  the  driving  signal  is  proportional  to  the  period  and  the 
averaging  volume  is  proportional  to  the  wavelength. 

Field  Methods 

The  pumping  system  uses  a  reservoir  on  the  surface,  a  pressurized  water  system, 
and  a  computer  controlled  servo-valve  to  supply  the  sinusoidally  varying  signal  (the  setup 
is  shown  in  Figure  5c).  We  chose  to  use  the  pneumatic  system  and  3  second  period  data 
for  the  geoprobe  source  because  of  logistics  and  time  critical  scheduling  of  the  geoprobe 
unit.  However,  it  seems  clear  that  the  pumping  system  could  have  been  adapted  equally 
well.  Using  two  vertical  multi-level  receiver  arrays  we  collected  data  simultaneously  in 
wells  HT-2  and  HT-3  while  using  the  Geoprobe  as  source.  The  field  setup  of  the 
Geoprobe  unit  is  shown  below  in  Figure  53a.  The  pneumatic  source  assembly  for  the 
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Geoprobe  unit  is  shown  in  Figure  53b  and  consisted  of  a  PVC  screened  one  foot  interval, 
blank  PVC  casing,  a  wafer  packer,  and  an  inflation  packer.  The  purpose  of  the  inflation 
packer  and  the  wafer  packer  was  to  isolate  the  Geoprobe  drill  string  from  the  source 
signal.  Unfortunately,  the  wafer  packer  burst  early  in  the  data  collection  so  only  the 
wafer  packer  was  available  for  isolation. 


Figure  53  a.  Geoprobe  unit  on  site  for  pneumatic  3  second  data  collection. 
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Figure  53b.  Pneumatic  source  assembly  consisting  of  a  PYC  screen  interval,  PYC  blank, 
wafer  packer,  and  inflation  packer. 

Data  Processing 

Data  processing  steps  for  these  data  were  performed  with  the  usual  series  of 
Visual  Basic  computer  programs  FitAmpPhase,  Fly draulicTom Anal,  and 
LeastSquaresSVD  discussed  and  used  earlier.  In  general,  analysis  for  this  study  used 
higher  quality  ZOP  data  to  estimate  K  at  relatively  discrete  locations  and  to  develop  a 
layered  model  grid  to  represent  the  aquifer  between  the  well  pairs.  Then  the  full  MOG 
data  sets  were  used  to  evaluate  aquifer  heterogeneity.  Linear  tomographic  inversion 
estimates  K,  which  is  derived  from  the  diffusivity  ratio  and  an  estimate  of  Ss.  Ss  is 
difficult  to  measure  and,  typically,  representative  values  are  usually  just  assumed  from 
literature  references  (Fetter,  2001).  A  corrected  Ss  estimated  from  baseline  HRST  K  and 
the  experimental  phase  data  were  used  during  the  modeling  and  inversion  process  to  best 
represent  the  aquifer  storage  characteristics  at  GEMS.  Some  records  exhibited  very  poor 
signal  to  noise  ratio  due  either  to  pressure  transducer  failure  or  attenuation  by  low  K 
layers.  In  these  cases  the  field  data  were  replaced  with  synthetic  data  or  surrogate  phase 
which  was  calculated  from  the  HRST  results  for  K.  Various  data  processing  techniques 
were  used  to  improve  the  fit  of  phase  data  to  the  layered  heterogeneous  models,  such  as 
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anisotropic  ratios  and  linear  constraint.  For  more  details  on  data  processing  see  Lyle 

(2011). 

Aquifer  heterogeneity  for  this  research  was  simulated  with  successive  model  runs 
using  ZOP  and  MOG  data  sets.  In  general,  each  of  the  well  pairs  had  three  modeling  and 
corresponding  inversion  steps  involving  horizontal  ZOP  data  to  develop  a  reduced  zone 
model  to  represent  the  aquifer,  MOG  data  to  evaluate  aquifer  anisotropy,  and  MOG  data 
to  evaluate  aquifer  lateral  heterogeneity.  First,  ZOP  data  were  modeled  on  a  mesh  grid 
spacing  that  corresponded  to  the  0.305  m  (1  ft)  spacing  of  the  source-receiver,  zero-offset 
profile  locations  within  the  well  pair.  Depending  on  the  number  of  source  intervals  in  a 
well  pair,  either  a  27  or  28  Zone  model  was  used  with  the  higher  quality  ZOP  data  to 
calculate  an  initial  K  for  each  source-receiver  location.  The  vertical  profde  of  this 
horizontal  K  data  was  used  to  develop  a  reduced  7  or  8  Zone  model  grid  which  better 
represents  the  expected  resolution  of  the  tomographic  method.  ZOP  ray  path  estimation 
through  the  7  or  8  Zone  model  and  corresponding  phase  inversion  determined  another 
vertical  distribution  of  K  values  which  were  subsequently  used  to  define  the  MOG  model 
runs.  Next,  the  reduced  7  or  8  Zone  model  along  with  ZOP  calculated  K  was  used  with 
the  MOG  data  to  simulate  isotropic  and  anisotropic  conditions. 

At  the  completion  of  the  successive  model  runs,  the  final,  constrained,  least 
squares  fit  K  values  were  contoured  against  elevation  and  radial  distance  between  source 
and  receiver  using  a  public  domain  program  called  QuickGrid.  The  program  contours 
between  points  written  in  an  x,y,z  format,  which  corresponds  to  radius,  elevation,  and  K 
values  determined  by  the  SVD  analysis  from  this  research.  In  all  of  the  contour  plots  of 
K,  the  source  well  is  on  the  left  side  and  receiver  well  is  on  the  right.  Using  the  HT-6  to 
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HT-3  well  pair  as  an  example  because  it  seemed  to  be  the  best  initial  data  set,  the 
modeling,  inversion,  statistical  K  determination,  anisotropy  evaluation  and  constraint 
factors  are  further  discussed  below.  For  a  more  complete  discussion  of  the  data 
processing  for  all  well  pairs  see  Lyle  (2011). 

ZOP  Modeling  and  Inversion 

To  obtain  the  expected,  optimal  1  m  (3.28  ft)  grid  resolution,  ZOP  data  were 
initially  modeled  on  a  finer  0.3048  m  (1  ft)  vertical  spacing  in  a  27  or  28  zone  model 
(Fig.  54).  This  zone  spacing  corresponds  to  the  experimental  source  and  receiver 
locations  in  a  well  and  the  number  of  zones  reflect  the  total  number  of  source  locations  in 
CPT.  ZOP  data  are  expected  to  have  the  best  quality  within  a  MOG  data  set  because  it 
has  the  shortest  ray  path  distance  between  the  source  and  receiver.  Consequently,  ZOP 
data  were  used  to  develop  a  coarser  model  grid  to  represent  the  aquifer  and  better  match 
the  expected  tomographic  resolution.  In  the  HT-6  to  HT-3  ZOP  model  below  (Fig.  54), 
source  locations  are  on  0.3048  m  (1ft)  centers  over  the  approximate  well  screen  interval 
of  232.0  to  240.2  m  msl  (761  to  788  ft  msl)  and  correspond  to  the  odd  numbered  nodes,  3 
-57,  within  the  model  grid.  Receiver  locations  are  directly  across  from  the  source  and 
correspond  to  the  even  number  nodes,  4  -  58.  Using  K,  Ss  and  ZOP  data, 
HydraulicTomAnal  modeled  the  28  horizontal  ray  paths,  each  through  their  respective 
zone,  to  generate  the  model  geometry  and  theoretical  model  phase  for  this  grid.  Either 
HRST  K  values  or  an  average  HRST  K  value  (0.003  ft/sec)  were  input  as  reasonable 
initial  constant  K  nodes  to  generate  the  model  matrix  and  phase. 
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Figure  54  -  Conceptual  28  zone  finite  difference  model  grid  for  HydraulicTomAnal. 


As  a  representative  data  set,  SVD  analysis  of  the  HT-6  to  HT-3  is  used  to  evaluate 
the  deterministic  K  calculated  from  the  ZOP  data  through  the  28  zone  model  (Fig.  54). 
The  K  values  correspond  to  the  28  zero-offset  source  and  receiver  locations  in  the  ZOP 
data  set  (Fig.  55).  The  profile  of  the  experimental  test  intervals  was  interpreted  into  a 
reduced  zone  model  with  thicker  vertical  zones  that  is  more  representative  of  the 
expected  tomographic  resolution.  Typically,  interpretation  of  the  ZOP  deterministic  K 
between  the  different  well  pairs  resulted  in  an  either  7  or  8  reduced-zone  model  (Fig.  56). 
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Figure  55  -  HT-6  to  HT-3  Deterministic  calculated  K  from  ZOP  data  set. 
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Figure  56.  HT-6  to  HT-3  interpretation  of  ZOP  deterministic  K. 
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The  thicker  vertical  elements  of  the  interpreted  8  Zone  model  (Fig.  56)  more 
closely  fall  within  the  expected  1  m  (3.28  ft)  resolution  of  the  tomographic  method.  The 
K  value  for  each  zone  is  an  average  of  the  deterministic  K  values  in  the  model  zone.  The 
HT-6  to  HT-3  example  model  has  eight  zones,  30  elements,  and  48  nodes  over  the 
approximate  232.0  to  240.2  m  msl  (761  to  788  ft  msl)  well  screen  interval  (Fig.  57). 

Each  zone  is  composed  of  either  two  elements  or  six  nodes  (Fig.  57).  The  center  nodes  in 
the  model  grid  define  two  lateral  elements;  these  extra  nodes  correspond  to  the  midpoint 
between  the  source  and  receiver  location.  These  center  nodes  will  be  useful  later  when 
we  consider  lateral  heterogeneity. 
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Figure  57.  Conceptual  8  zone  model  grid  with  8  zones  (blue)  and  30  elements  (red). 

To  check  if  the  7  or  8  zone  model  adequately  represents  the  larger  ZOP  data  set, 
the  modeled  reduced  zone  model  phase  was  plotted  against  the  ZOP  experimental  phase 
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(Fig.  58).  At  HT-6  to  HT-3,  the  smoother  curve  of  the  8  zone  model  phase  reflects  the 
average  value  of  the  deterministic  K  assigned  to  the  thicker  layers  of  the  model.  The 
sharper  curve  of  the  experimental  phase  reflects  the  experimental  phase  at  each  of  the 
initial  layers  used  in  the  28  zone  model. 


HT-6  to  HT-3  Node  Data  Summary 
8  Zone  K/ve  Model  Phase  vs.  Exp.  Phase 


Figure  58.  Experimental  ZOP  phase  shift  plotted  against  ZOP  8  zone  model  phase  shift. 

SVD  analysis  of  the  reduced  zone  model  with  the  ratio  of  zones  to  rays  no  longer 
1 : 1  generates  a  least  squares  fit  of  K.  The  reduced  zone  calculated  K  is  an  average  value 
of  the  deterministic  K  in  each  zone  (Fig.  59).  Up  to  this  point  SYD  analysis  and 
modeling  used  an  Ss  estimate  of  0.00001  to  derive  K  from  diffusivity.  The  initial  IE-05 
value  is  considered  to  be  representative  of  the  aquifer  storage  characteristic  at  GEMS 
based  on  literature  references  and  previous  work  by  Wachter  (2008)  and  others.  The  Ss 
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value  and  K  from  this  inversion  was  further  corrected  to  more  closely  match  the  baseline 
HRST  K  data  and  the  phase  shift  data  before  further  modeling  with  the  MOG  data  sets. 
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Figure  59.  28  Zone  deterministic  K  plotted  against  8  Zone  calculated  K. 


Ss  and  K  were  corrected  to  match  GEMS  site  conditions.  Because  phase  should 

vary  linearly  with  ,  the  corrected  Ss  value,  S scorr,  was  calculated  from  the  average 

HSRT  K  values,  A^HRSTAve,  and  the  average  deterministic  K  values  from  the  28  zone  ZOP 
model,  KDet  Zop_Ave  ■  Based  on  the  theory,  all  analyses  were  simply  adjusted  to  the  new 
Ss corr  value  with  a  correction  factor  derived  from  a  ratio  of  the  average  slug  test  and 
model  deterministic  K  values: 
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(12) 


^ HRST  Ave  ^  ,•  ^  , 

- = - =  Correction  Factor 

K-Det_ZOP_Ave 

When  phase  shift  is  constant,  the  relationship  between  Ss  and  K  in  the  basic  phase  shift 
equation  is  linear,  so  changes  to  K  or  Ss  will  vary  in  proportion  to  each  other. 
Accordingly,  the  deterministic  and  least  squares  calculated  K  values  were  multiplied  by 
the  correction  factor  to  reflect  the  corrected  Ss.  The  corrected  HT-6  to 
HT-3  deterministic  K  values  from  the  reduced  zone  model  are  presented  below  (Fig.  60). 
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Figure  60.  HT-6  to  HT-3  correction  for  Ss  and  K. 


A  comparison  of  the  corrected  vertical  K  profile,  which  best  represents  the  storage 
characteristic  of  the  aquifer,  to  the  uncorrected  values  is  presented  in  Figure  61. 
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Figure  61.  HT-6  to  HT-3  vertical  K  profile  adjusted  for  the  calculated  Ss  that  represents 
site-specific  aquifer  storage  characteristics. 


MOG  Modeling,  Inversion,  and  Anisotropy 

After  the  7  or  8  reduced-zone  model  geometry  was  developed  and  verified,  and 
the  Ss  and  calculated  K  were  corrected  to  site-specific  conditions,  modeling  with  MOG 
data  sets  was  initiated  with  the  reduced-zone  model  to  determine  a  representative  ratio  of 
anisotropy  for  data  processing.  In  an  idealized  system,  isotropic  conditions  are  often 
assumed  where  Kh  =  Ky.  In  a  natural  system  anisotropy  prevails  and  typically  Kh  »> 
Ky  (Domenico  and  Schwartz,  1998).  The  anisotropic  ratios  applied  to  the  model  simulate 
a  more  realistic  flow  regime  to  better  imitate  the  directional  dependency  exhibited  by 
alluvial  sediments  oriented  by  grain  size  and  direction.  In  this  case,  an  anisotropic  ratio 
of  10  is  KH  =  10KV  and  an  anisotropic  ratio  of  2  is  KH  =  2KV.  Anisotropic  evaluation 
was  completed  for  both  the  pumped  hydraulic  CPT  data  (30-sec  period)  and  the 
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pneumatic  CPT  data  (3-sec  period).  Different  anisotropic  ratios  were  evaluated  with 
HydraulicTomAnalV21Aniso  which  functions  the  same  as  HydraulicTomAnalV21,  but 
generates  anisotropic  effects  across  an  element  or  zone.  Anisotropic  ratios  can  be  applied 
in  multiple  combinations  over  different  layers  to  simulate  greater  or  lesser  anisotropy.  A 
typical,  complete  CPT  well  pair  may  have  28  MOG  data  sets  with  784  rays  instead  of  just 
the  28  ZOP  rays  initially  used  to  develop  the  reduced-zone  model  of  the  aquifer  for 
heterogeneity  evaluation.  With  the  same  K  inputs  between  the  anisotropic  and  isotropic 
models,  varying  anisotropy  can  be  statistically  evaluated  after  model  inversion.  The 
offset  MOG  rays  should  measure  the  anisotropic  variation  in  the  aquifer,  so  an 
anisotropic  correction  should  theoretically  improve  the  data  fit  to  the  model. 

Considerable  experimentation  with  anisotropy  ratios  was  carried  out  (Lyle,  2011). 
It  was  anticipated  that  a  single  best  anisotropic  ratio  could  be  found  and  applied  to  all 
well  pairs.  However,  a  significant  statistical  deviation  occurred  when  processing  the  3- 
sec  data  sets,  so  additional  anisotropic  evaluation  was  completed  for  the  pneumatic  3-sec 
data.  Evaluation  results  suggested  that  different  anisotropic  ratios  and  data  constraints 
are  needed  to  adequately  model  CPT  data  with  different  oscillating  periods.  So,  after  the 
initial  assessment,  different  anisotropy  ratios  were  applied  to  the  30  and  3-sec  data  sets, 
but  all  the  model  layers  used  a  single,  best  case  anisotropy  ratio  (e.g.,  2  or  10). 

Inversion  and  least  squares  fitting  to  the  reduced  zone  model  was  evaluated  under 
different  constraints  (weights  on  initial  estimates  of  K),  as  well  as  isotropic  and 
anisotropic  ratios  of  2  and  10  for  the  two  CPT  sources.  Chi  squared  and  standard 
deviation  values  from  the  least  squares  fit  for  the  hydraulic  30-sec  data  indicate  that  an 
anisotropic  ratio  of  10  produced  a  poorer  data  fit  (representative  data  for  HT-1  to  HT-3 
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shown  in  Fig.  62).  Both  the  isotropic  and  the  anisotropic  2  models  were  relatively  good, 
although  the  anisotropic  model  provided  a  slightly  better  data  fit.  Most  of  the  K  values 
solved  from  the  anisotropic  2  model  were  near  or  less  than  10%  standard  deviation. 


HT-1  to  HT-3 

Least  Squares  Fit  K  -  %  Stdr  Dev  K 
8  Zone  Iso,  Aniso  2,  &  Aniso  10 


Figure  62.  HT-1  to  HT-3  least  squares  fit  K  under  isotropic,  anisotropic-2,  and 
anisotropic- 10  conditions. 


The  overlap  of  field  data  and  the  calculated  phase  values  (e.g.,  isotropic, 
anisotropic  2,  anisotropic  10,  etc.)  should  indicate  the  relative  goodness  of  fit  between  the 
experimental  data  and  the  straight  ray  approximation.  As  an  example,  MOG  data  from 
the  first  CPT  source  location  at  HT-1  to  HT-3  for  the  isotropic,  anisotropic  2,  and 
anisotropic  10  scenarios  are  presented  below  (Fig.  63  to  65).  In  this  example,  the  CPT 
source  location  is  near  the  bottom  of  HT-3  on  the  left  and  the  28  receiver  locations  in 
HT-1  are  on  the  right.  The  ray  path  of  the  CPT  must  travel  farther  to  reach  the  receivers 
located  in  the  higher  parts  of  the  well  screen,  so  the  phase  shift  will  increase  with 
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distance  and  form  a  half  parabola  at  this  location.  The  curves  most  closely  match  with  an 
anisotropy  ratio  of  2  (Fig.  64),  which  correlates  with  the  data  fit  indicated  by  the  percent 
standard  deviation  results  displayed  in  Figure  62.  Additional  experiments  were  run 
testing  different  anisotropy  ratios  for  high  K  and  low  K  zones,  but  no  significant 
improvement  was  noted.  So,  an  anisotropic  ratio  of  2  was  chosen  to  evaluate  all  the 
MOG  data  sets  for  the  30-sec  CPT  data  (i.e.,  the  pumped  hydraulic  source). 


HT-1  to  HT-3 
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Figure  63.  HT-1  to  HT-3  experimental  phase  shift  and  SVD  calculated  phase  shift  under 
isotropic  conditions. 
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Figure  64.  HT-1  to  HT-3  experimental  phase  shift  and  SVD  calculated  phase  shift  with 
an  anisotropic  ratio  2. 
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Figure  65.  HT-1  to  HT-3  experimental  phase  shift  and  SVD  calculated  phase  shift  with 
an  anisotropic  ratio  10. 
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Conversely,  modeling  for  the  3-sec  pneumatic  data  from  HT-GP  to  HT-2  using  an 
anisotropic  ratio  of  2  did  not  result  in  the  best  data  fit.  Large  chi  squared  and  standard 
deviation  values  resulted  from  the  inversion  of  this  anisotropic  scenario  (Fig.  66).  In 
addition,  some  of  the  fitted  K  values  were  out  of  the  expected  range,  suggesting  some 
constraint  may  be  needed.  Inverse  problems  are  sometimes  constrained  by  other  sources 
of  data  or  by  mathematical  methods  if  the  results  prove  to  be  unrealistic.  In  some  of  the 
earlier  tomographic  research  on  this  project,  Wachter  (2008)  found  that  some  of  the 
calculated  K  values  were  an  order  of  magnitude  higher  than  the  rest  of  the  data  set  and  a 
weighting  factor  within  LeastSquaresSVD  program  was  used  to  constrain  the  inversion 
closer  to  site-specific  HRST  K  values.  Up  to  this  point  in  the  data  processing,  all  the 
ZOP  and  anisotropy  model  inversions  were  completed  with  a  constraint  factor  (CLS)  of 
0,  which  is  unconstrained.  Increasing  constraint  values  gives  more  weight  to  input  K 
values  and  therefore  restrict  deviations  of  the  inversion  results  from  the  input  model  K 
values. 
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Figure  66.  HT-GP  to  HT-2  least  squares  fit  K  under  isotropic,  anisotropic-2,  and 
anisotropic- 10  CLS  1  conditions. 


The  3 -sec  data  from  the  pneumatic  HT-GP  to  HT-2  well  pair  were  remodeled 
with  isotropic,  anisotropic  2,  and  anisotropic  10  models  to  determine  a  better  aquifer 
model  and  inversion  for  the  pneumatic  data.  A  constraint  factor  of  1  on  the  anisotropic 
10  model  was  chosen  for  SVD  analysis,  which  still  lends  equal  weight  to  both  the  input  K 
of  the  model  and  the  inverted  results  to  help  avoid  an  artificial  data  fit  to  the  model.  The 
data  fit  was  much  improved;  the  constraint  factor  of  1  kept  the  K  values  in  the  expected 
range. 


Interestingly,  the  larger  anisotropy  ratio  improved  the  overall  data  fit  of  the  offset 
diagonal  rays  in  the  MOG  pneumatic  data.  Initial  inversion  of  the  anisotropic  2  model 
resulted  in  a  relatively  good  fit  between  the  experimental  phase  shift  data  and  the  model 
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generated  phase  shift  for  many  of  the  MOGs  at  HT-GP  to  HT-2,  although  the  fit  of  the 
offset  rays  increasingly  deteriorated  with  distance  (Fig.  67).  This  was  somewhat 
expected  for  the  short  period  data  sets  since  higher  frequencies  are  more  strongly 
attenuated.  Attenuation  of  the  signal  is  inversely  proportional  to  K  of  the  medium  and 
the  period,  so  the  factor  of  10  difference  between  the  two  CPT  the  periods  and  low  K 
alluvium  in  the  upper  portion  of  the  test  interval  will  exacerbate  the  attenuation.  It  was 
unknown  if  the  attenuation  with  distance  through  low  K  material  would  negate  the 
greater  resolution  gained  by  the  short  period  CPT  source,  but  data  analysis  of  MOG 
phase  data  at  HT-GP  to  HT-2  suggests  that  a  greater  anisotropy  ratio  seems  to  help  fit  the 
low  K  data  points  (i.e.,  high  phase)  of  the  long  offset  rays.  The  calculated  anisotropic  10 
phase  from  the  longest  rays  or  uppermost  receivers  more  closely  fits  the  experimental 
data  set  than  the  calculated  anisotropic  2  phase  (Fig.  68  vs  Fig.  67).  The  30-sec  period 
data  sets  were  relatively  insensitive  to  different  anisotropy  ratios  greater  than  2,  so  this 
suggests  that  the  3 -sec  CPT  period  data  are  more  sensitive  to  aquifer  anisotropy  than  the 
30-sec  CPT  period  data  in  spite  of  signal  attenuation  associated  with  increasing  radial 
distance,  decreasing  K,  and  decreasing  period. 
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Figure  67.  Experimental  field  phase  shift  and  phase  shift  calculated  using  unconstrained 
SVD  for  the  first  MOG  from  HT-GP  to  HT-2. 
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Figure  68.  Experimental  field  phase  shift  and  phase  shift  calculated  using  constrained 
SVD  for  the  first  MOG  from  HT-GP  to  HT-2. 
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Most  of  the  zones  from  the  constrained  SVD  analysis  of  the  anisotropic  10  model 
at  HT-GP  to  HT-2  had  an  approximately  10%  standard  deviation  value,  although  two 
zones  in  the  basal,  higher  K  portion  of  the  aquifer  were  somewhat  higher,  but  still 
remained  less  than  20%  (Fig.  66).  Based  on  the  statistical  evaluation  of  the  SVD 
inversion,  an  anisotropic  ratio  of  10  and  constrained  least  squares  factor  of  1  was  chosen 
to  evaluate  all  the  MOG  data  sets  with  a  3-sec  oscillation  period  (i.e.,  pneumatic  source) 
for  lateral  heterogeneity.  Lateral  heterogeneity  modeling  and  inversion  constraint  are 
further  discussed  in  the  subsequent  section. 

MOG  Modeling  with  Lateral  Heterogeneity 

Finally,  the  effect  of  lateral  heterogeneity  was  evaluated  with  a  24  zone  model 
(Fig.  69).  Each  of  the  horizontal  zones  that  correspond  to  the  interpreted 
hydrostratigraphic  layers  of  the  GEMS  aquifer  (Fig.  56)  was  subdivided  to  include  three 
lateral  zones.  The  HT-6  to  HT-3  example  model  has  24  zones,  30  elements,  and  48 
nodes  over  the  approximate  232.0  to  240.2  m  msl  (761  to  788  ft  msl)  well  screen  interval. 
The  nodes  at  the  source  well  (e.g.,  node  1  and  4)  and  the  midpoint  between  the  well  pair 
(e.g.,  node  2  and  5)  and  receiver  well  (node  3  and  6)  define  the  lateral  zones.  Instead  of 
holding  K  constant  across  a  single  horizontal  zone,  the  K  value  for  each  layer  was 
allowed  to  vary  laterally  at  these  nodal  points.  HydraulicTomAnalV21  Aniso  linearly 
interpolates  between  the  nodes  to  simulate  the  effects  of  lateral  heterogeneity  across  the 
model  (Fig.  69).  As  determined  by  the  MOG  anisotropy  evaluation,  phase  data  were 
simulated  through  the  aquifer  with  an  anisotropy  ratio  of  2  for  the  30-sec  data  and  an 
anisotropy  ratio  of  10  for  the  3-sec  data. 
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Figure  69.  Conceptual  24  zone  model  grid  for  determining  lateral  heterogeneity. 


The  inversion  and  least  squares  fit  to  the  lateral  heterogeneity  model  for  the  30- 
sec  and  3-sec  period  CPT  data  was  evaluated  under  different  constraints.  As  discussed  in 
the  previous  section,  MOG  ray  path  estimation  was  generated  through  the  lateral 
heterogeneity  model  for  the  30  and  3-sec  data  sets  with  anisotropy  ratios  of  2  and  10, 
respectively.  The  lateral  subdivision  of  the  reduced  zone  model  adds  more  variables  to 
the  SVD  inversion,  possibly  increasing  the  generation  of  non-unique  results  or 
uncertainty  within  the  data  fit.  SVD  constraint  factor  experiments  for  the  30  and  3 -sec 
CPT  data  sets  were  performed. 

Chi  squared  and  standard  deviation  values  for  the  calculated  K  from  the  30-sec 
data  indicate  that  a  constraint  factor  of  1  is  generally  acceptable  for  the  inversion  of  the 
pumped  hydraulic  CPT  data  (Fig.  44  and  Fig.  45).  Most  of  the  standard  deviation  values 
are  approximately  10%.  However,  the  fit  is  less  good  in  the  basal,  high  K  portions  of  the 
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aquifer  although  the  percent  standard  deviation  does  remain  below  20%.  The  calculated 
K  values  from  the  30-sec  MOGs  are  still  within  the  range  typically  encountered  at 
GEMS,  so  a  constraint  factor  of  1  was  chosen  to  evaluate  lateral  heterogeneity  for  the  30- 
sec  CPT  data. 


HT-1  to  HT-3 

Least  Squares  Fit  K  -  %  Stdr  Dev  K 
24  Zone  Ansio  2  -  CLS  1  vs  CLS  10 


Figure  70.  HT-1  to  HT-3  least  squares  fit  K  under  anisotropy-2  CLS  1  (Ave  10.73  %) 
and  anisotropy-2  CLS  10  (Ave  5.72%)  conditions. 
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HT-6  to  HT -3 

Least  Squares  Fit  K  -  %  Stdr  Dev  K 
24  Zone  Ansio  2  -  CLS  1  vs  CLS  10 


Figure  71.  HT-6  to  HT-3  least  squares  fit  K  under  anisotropy-2  CLS  1  (Ave  9.95%)  and 
anisotropy-2  CLS  10  (Ave  5.23%)  conditions. 

Initially,  a  constraint  factor  of  1  was  used  for  the  3 -sec  CPT  data  inversion. 
However,  once  the  HT-GP  to  HT-2  data  were  processed,  the  chi  squared  and  percent 
standard  deviation  values  of  the  calculated  K  from  the  SVD  analysis  increased  markedly 
(33.6  %  vs  the  approximate  10%  average  of  the  hydraulic  data),  and  indicated  that  the 
pneumatic  data  were  responding  differently  to  the  inversion  and  fit  to  the  lateral 
heterogeneity  model  (Fig.  2).  An  additional  SVD  inversion  with  a  constraint  factor  of  10 
was  completed  to  evaluate  the  pneumatic  data  at  HT-GP  to  HT-2.  Some  error  was 
removed  with  the  additional  constraint  and  the  average  percent  standard  deviation  was 
reduced  to  about  12%,  although  a  plot  of  the  percent  standard  deviation  has  little 
variability,  so  the  data  fit  appears  somewhat  artificially  constrained.  The  calculated  K 
values  with  a  least  squares  constraint  of  1  are  still  within  the  expected  range  of  reported 
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K  values  at  GEMS  so  this  constraint  factor  was  still  assumed  for  this  CPT  well  pair. 
Again,  generally,  the  fit  is  less  good  in  the  basal  portion  of  the  aquifer. 

In  contrast  to  the  other  hydraulic  and  pneumatic  CPT  well  pairs,  the  remaining 
pneumatic  location,  HT-GP  to  HT-3,  required  a  constraint  factor  of  10  to  generate  K 
values  within  the  expected  range  at  GEMS  (Fig.  73).  Although  a  constraint  factor  of  1 
still  resulted  in  a  slightly  better  data  fit  than  the  other  pneumatic  well  pair,  (29.2%  vs 
33.6%  average  standard  deviation)  the  K  values  deviated  from  the  expected  range  and,  in 
particular,  one  K  value  from  Zone  19  (0.0624  ft/sec)  was  an  order  of  magnitude  greater 
than  the  rest  of  the  data  set  with  a  percent  standard  deviation  error  of  243%  (Fig.  73).  A 
constraint  factor  of  10  resulted  in  reasonable  K  values  and  percent  standard  deviation 
(16.9%);  therefore,  it  was  used  for  lateral  heterogeneity  evaluation  at  HT-GP  to  HT-3. 


HT-GP  to  HT-2 

Least  Squares  Fit  K  -  %  Stdr  Dev  K 
21  Zone  Aniso  2  &  10  -  CLS  1  vs  CLS  10 


Figure  72.  HT-GP  to  HT-2  least  squares  fit  K  under  anisotropic-2  CLS  1  (Ave  33.6%), 
anisotropic- 10  CLS  1  (Ave  26.5%),  anisotropic- 10  CLS  10  (Ave  12.1%)  conditions. 
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HT-GP  to  HT-3 

Least  Squares  Fit  -  %  Stdr  Dev  K 
24  Zone  Aniso  2  &  10  -  CLS  1  vs  CLS  10 


Figure  73.  HT-GP  to  HT-3  least  squares  fit  K  under  anisotropic-2  CLS  1  (Ave  45.5%), 
anisotropic- 10  CLS  1  (Ave  29.2%),  anisotropic- 10  CLS  10  (Ave  16.9%)  conditions. 


Hydraulic  Conductivity  Distributions  -  Pumped  Hydraulic  CPT  Source 

After  the  data  processing  outlined  above  was  completed  for  all  well  pairs,  contour 
plots  were  made  of  K  values  plotted  against  elevation  and  the  radial  distance  between  a 
well  pair  using  the  graphing  program,  QuickGrid.  The  program  contours  between  points 
written  in  an  x,y,z  format,  which  corresponds  to  the  radius,  elevation,  and  constrained  K 
value  determined  from  SVD  analysis.  Five  pumped  hydraulic  CPT  well  pairs  were 
plotted  with  the  receiver  well  on  the  left  (HT-3)  and  the  source  wells  on  the  right  (HT-1, 
HT-2,  HT-4,  HT-5  and  HT-6).  The  K  values  from  this  radial  array  were  obtained  with 
the  pumped  hydraulic  CPT  source  which  had  a  30-sec  oscillating  period.  A  best  case 
anisotropic  ratio  and  inversion  constraint  factor  for  the  lateral  model  grid  and  SVD 
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inversion,  respectively,  were  chosen  based  on  the  evaluation  of  different  anisotropic 
scenarios  of  the  reduced  zone  MOG  model  at  HT-3  to  HT-6  and  HT-3  to  HT-1,  as  well 
as,  the  expected  range  of  K  values  at  GEMS  (0.0003  to  0.003  m/sec  [0.001  to  0.01 
ft/sec]).  In  particular,  the  well  pair  at  HT-3  to  HT-6  was  thought  to  have  the  best  overall 
data  quality  and  was  used  as  a  benchmark  reference  for  the  other  data  sets.  K  evaluation 
by  individual  well  pair  is  further  discussed  below.  Tables  of  K  values  are  available  in 
Lyle  (2011). 

The  K  values  for  the  CPT  well  pair  HT-3  to  HT-6  are  presented  in  Fig.  74  and  are 
within  the  range  expected  at  GEMS  (0.0003  to  0.0018  m/sec  [0.0010  to  0.0058  ft/sec]). 
The  contour  trend  follows  the  expected  results  for  the  GEMS  lithology  and  HRST  results, 
with  high  K  values  in  the  coarser,  basal  portion  of  the  aquifer  and  low  K  values  in  the 
finer,  upper  portion  of  the  aquifer.  The  HT-3  to  HT-6  data  set  did  not  have  any  poor 
quality  data  and  the  data  fit  (9.95%  average  standard  deviation  K)  was  relatively  good. 
Some  of  the  larger  percent  standard  deviation  occurred  in  the  high  K  portions  of  the  test 
interval. 

The  K  values  presented  in  Figure  75  for  the  well  pair  HT-3  to  HT-1  (0.0002  to 
0.0017  m/sec  [0.0008  to  0.0056  ft/sec])  are  generally  within  the  range  expected  at 
GEMS,  although  the  low  range  K  value  is  slightly  lower  than  what  has  been  reported  in 
the  past  for  the  fine  grained  portion  of  the  GEMS  aquifer.  However,  the  slight  deviation 
is  not  considered  significant  and  the  overall  contour  trend  of  the  plot  still  follows  the 
expected  GEMS  lithology  and  HRST  results.  The  HT-3  to  HT-1  data  set  did  not  have 
any  HRST  replacement  data  used  and  the  data  fit  (10.7%  average  standard  deviation  K) 
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was  relatively  good.  As  usual,  some  of  the  larger  percent  standard  deviation  occurred  in 
the  high  K  portions  of  the  test  interval. 

The  data  set  presented  in  Figure  76  from  HT-3  to  HT-2  has  a  range  of  K  values 
(0.0002  to  0.0016  m/sec  [0.0005  to  0.0051  ft/sec])  that  is  generally  within  the  expected 
range  at  GEMS  and  the  contour  plot  depicts  the  expected  K  distribution.  The  well  pair 
did  not  have  any  HRST  replacement  data  but  the  data  set  was  not  quite  as  good  as  the 
benchmark  well  pairs.  The  amount  of  error  between  calculated  and  observed  phases 
(14.6%  average  standard  deviation  K)  was  greater  than  the  benchmark  well  pairs  at  HT-6 
to  HT-3  and  HT-1  to  HT-3.  Some  of  the  error  could  likely  be  removed  from  the  data  set 
with  greater  constraint  during  SVD  analysis,  but  because  the  K  trends  remain  reasonable, 
additional  constraint  was  deemed  unnecessary. 

The  CPT  data  presented  in  Figure  77  from  HT-3  to  HT-4  include  some  HRST 
replacement  data  in  the  middle  of  the  test  section  (see  Lyle,  2011  for  details).  The  ranges 
of  K  values  (0.0002  to  0.0019  m/sec  [0.0006  to  0.0061  ft/sec])  are  within  the  expected 
range  and  the  contour  plot  depicts  a  reasonable  K  distribution,  although  some  of  the 
heterogeneity  graduation  seen  between  the  lower  and  upper  zones  in  the  benchmark  well 
pairs  seems  to  be  suppressed  or  absent  in  the  zones  with  replacement  data.  The  K  values 
are  still  within  the  range  of  other  K  data  and  the  amount  of  error  between  calculated  and 
observed  phases  is  relatively  good  (9.43%  average  standard  deviation),  so  the  plot 
appears  to  be  a  representative  depiction  of  the  aquifer  heterogeneity  expected  at  GEMS. 
Comparison  of  this  plot  to  Wachter’s  (2008)  interpretation  of  the  well  pair  would  tend  to 
confirm  this  lack  of  heterogeneity  in  this  portion  of  the  aquifer.  The  upper  and  middle 
zones  are  largely  low  K  alluvium  with  the  high  K  alluvium  limited  to  just  the  basal 
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portion  of  the  aquifer,  although  the  K  values  in  the  basal  portion  from  this  research  are 
somewhat  lower  in  magnitude  and  not  as  consistently  extensive  across  the  section. 

The  CPT  data  shown  in  Figure  78  from  HT-3  to  HT-5  include  replacement  HRST 
data  for  over  half  of  the  upper  test  interval  (see  Lyle,  2011  for  details),  due  to  pressure 
transducer  failures.  The  ranges  of  K  values  (0.0002  to  0.0024  m/sec  [0.0005  to  0.0079 
ft/sec])  are  within  the  expected  range  at  GEMS.  Compared  to  the  other  hydraulic  CPT 
well  pair  with  HRST  replacement  data  (HT-4  to  HT-3),  the  data  fit  is  not  as  good  (1 1.6% 
vs  9.43%  average  standard  deviation).  However,  the  amount  of  error  between  calculated 
and  observed  phases  at  this  well  pair  (11 .6%  average  standard  deviation  K)  is  still 
consistent  with  HT-1  to  HT-3  (10.7%  average  standard  deviation  K),  which  is  one  of  the 
benchmark  well  pairs.  The  greatest  standard  deviation  again  appears  in  the  lower  high  K 
zone.  The  heterogeneity  trends  are  consistent  with  the  other  well  pairs,  but  a  higher  K 
zone  in  the  upper  portion  of  the  aquifer  near  the  source  well  (HT-3)  seems  somewhat 
more  pronounced  and  probably  reflects  the  direct  measurement  of  the  HRST  K  point 
source  data  used  to  replace  the  bad  data  for  the  CPT  well  pair. 

The  pumped  hydraulic  CPT  data  (30-sec  period)  collected  from  the  radial  well 
area  replicate  the  aquifer  interval  that  Wachter  (2008)  tested  with  a  pneumatic  CPT 
source  (3-4  second  period).  The  summary  results  of  that  earlier  research  were  evaluated 
and  compared  to  the  present  research  to  evaluate  the  difference  between  the  two  different 
CPT  source  methods.  Wachter,  to  identify  the  best  model  to  evaluate  heterogeneity  as 
well  as  inversion  constraint,  used  a  number  of  different  model  configurations,  ray  paths 
and  different  Ss  values  (IE-05  and  1.5E-05).  The  phase  depends  on  a  ratio  between  Ss 
and  K,  so  changes  in  Ss  will  also  result  in  changes  in  K.  This  introduces  a  potential 
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source  of  error,  due  to  the  difficulty  of  measuring  Ss  in  situ.  In  this  research,  an  attempt 
to  limit  this  error  was  implemented  by  correcting  Ss  using  the  measured  phase  and  HRST 
or  ZOP  K.  The  corrected  Ss  was  used  for  modeling  and  inversion.  Wachter  determined 
that  a  Ss  value  of  1.5E-05  and  about  750  rays  produced  the  best  results.  The  corrected  Ss 
values  (about  1.1E-05  to  1.69E-05)  and  MOG  rays  (756  to  784)  used  for  this  research 
were  comparable  to  Watcher’s  earlier  tomographic  work.  Wachter’ s  model 
configurations  were  slightly  different,  which  used  elements  instead  of  nodes  for  straight 
ray  approximation  through  the  model  grid,  and  his  model  generally  had  fewer  total  zones 
(e.g.,  16  vs  24).  The  different  model  configuration  for  the  present  research  was 
developed  from  the  use  of  the  relatively  unattenuated  ZOP  data  to  create  a  representative 
aquifer  model.  It  was  expected  that  utilizing  a  node  model  grid  along  with  more  vertical 
and  lateral  zones  and  initial  use  of  ZOP  phase  data  could  improve  the  resolution  of  the 
method. 

The  range  of  K  values  derived  by  Wachter  (2008)  with  the  4-second  CPT  period 
are  relatively  consistent  with  the  ones  in  this  research  and  follow  the  expected  trends 
based  on  the  aquifer  lithology  and  range  of  reported  K  values  at  GEMS.  However,  in 
general,  there  was  some  difference  in  heterogeneity  resolution  between  the  two  CPT 
sources.  Wachter,  with  the  exception  of  the  HT-4  to  HT-3  well  pair,  was  able  to  resolve 
a  somewhat  higher  K  zone  (e.g.,  0.0015  -  0.0045  m/sec  [0.0049  -  0.0148  ft/sec])  in  the 
middle  of  the  aquifer.  In  contrast,  the  K  values  derived  with  the  pumped  hydraulic  CPT 
were  mostly  lower  in  the  intermediate  zone  (e.g.,  0.0006  to  0.001  m/sec  [0.0020  to 
0.0033  ft/sec]).  Although,  as  an  exception,  the  HT-5  to  HT-3  well  pair  did  have 
somewhat  higher  K  values  (e.g.,  0.0015  m/sec  [0.0049  ft/sec])  through  this  zone. 
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This  heterogeneity  difference  through  the  middle  zone  could  reflect  the  way  the 
data  plots  were  constrained  during  contouring  or  a  difference  between  the  periods  of  the 
two  CPT  sources.  Wachter’s  contour  plots  used  SVD-determined  K  values  which  were 
constrained  by  HRST  K  data  along  the  source  and  receiver  well  locations.  The  contour 
plots  for  this  research  were  not  constrained  by  HRST  K.  The  HT-5  to  HT-3  contour  (Fig. 
58)  may  reflect  this  constraint  difference.  This  well  pair  had  a  significant  portion  of  the 
field  data  replaced  with  surrogate  phase  data  which  are  estimated  from  HRST  K,  so  the 
SVD  analysis  there  may  be  replicating  some  of  the  HRST  K  constraint  within  the  3-4 
second  CPT  period  data  sets.  Also,  in  consideration  of  the  period  differences,  results 
from  the  numerical  modeling  of  the  heterogeneity  extension  indicate  that  the  30-second 
CPT  period  does  not  have  as  much  resolution  as  the  3-4  second  CPT  period,  so  some  of 
this  heterogeneity  loss  is  possibly  due  to  the  different  CPT  source  periods. 
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Figure  74.  K  values  from  constrained  SVD  analysis  of  784  rays  at  HT-3  (receiver)  to 
HT-6  (source).  K  values  range  from  0.0003  to  0.0018  m/sec  (0.0010  to  0.0058  ft/sec). 
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Figure  75.  K  values  from  constrained  SVD  analysis  of  756  rays  at  HT-3  (receiver)  to 
HT-1  (source).  K  values  range  from  0.0002  to  0.0016  m/sec  (0.0005  to  0.0051  ft/sec). 
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Figure  76.  K  values  from  constrained  SVD  analysis  of  756  rays  at  HT-3  (receiver)  to 
HT-2  (source).  K  values  range  from  0.0002  to  0.0016  m/sec  (0.0005  to  0.0051  ft/sec). 
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Figure  77.  K  values  from  constrained  SVD  analysis  of  784  rays  at  HT-3  (receiver)  to 
HT-4  (source).  K  values  range  from  0.0002  to  0.0019  m/sec  (0.0006  to  0.0061  ft/sec). 
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Figure  78.  K  values  from  constrained  SVD  analysis  of  784  rays  at  HT-3  (receiver)  to 
HT-5  (source).  K  values  range  from  0.000  to  0.00  m/sec  (0.000  to  0.00  ft/sec). 


140 


Hydraulic  Conductivity  Distributions  -  Pneumatic  CPT  Geoprobe  Source 

Contour  plots  were  made  of  K  values  plotted  against  elevation  and  the  radial 
distance  between  the  Geoprobe  source  location  and  two  receiver  wells  using  the  graphing 
program,  QuickGrid.  Two  pairs  of  results  for  K  were  plotted  with  the  source  location 
(HT-GP)  on  the  left  and  receiver  wells  on  the  right  (i.e.,  HT-2  or  HT-3).  The  K  values 
from  this  triangular  array  were  obtained  with  the  pneumatic  CPT  method  at  the  source 
and  a  3 -sec  oscillating  period.  A  best  case  anisotropic  ratio  and  inversion  constraint 
factor  for  the  lateral  model  grid  and  SVD  inversion,  respectively,  were  chosen  based  on 
the  anisotropic  evaluation  of  the  reduced-zone  MOG  model  and  the  lateral  heterogeneity 
model  at  HT-GP  to  HT-2.  Data  fit  to  these  models  was  significantly  improved  with  a 
constraint  factor  of  1  (which  better  fit  the  replacement  data)  and,  by  increasing  the 
anisotropic  ratio  from  2  to  10  (which  better  fit  the  farthest  offset  rays  of  the  MOG).  This 
improved  fit  due  to  the  application  of  a  larger  anisotropic  ratio  was  not  apparent  in  the 
30-sec  CPT  data  sets  and  suggests  that  the  3-sec  data  sets  were  more  sensitive  to 
anisotropy.  However,  as  discussed  previously,  additional  inversion  constraint  was 
needed  for  the  HT-GP  to  HT-3  well  pair  to  suppress  a  K  data  point  that  was  an  order  of 
magnitude  greater  than  the  rest  of  the  K  values  in  the  lateral  model.  These  constraints, 
anisotropy,  and  the  expected  range  of  K  values  at  GEMS  (0.0003  to  0.003  m/sec  [0.001 
to  0.0098  ft/sec])  were  used  to  evaluate  the  lateral  heterogeneity.  K  evaluation  by 
individual  well  pair  is  further  discussed  and  presented  below. 

The  K  values  presented  in  Figure  79  from  HT-GP  to  HT-2  involve  some 
replacement  data  in  the  middle  and  upper  portion  of  the  test  section,  replacing  some 
unusable  field  data  (see  Lyle,  201 1  for  details).  The  ranges  of  K  values  (0.0003  to  0.002 
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m/sec  [0.0009  to  0.0067ft/sec])  are  within  the  expected  range  at  GEMS.  Different 
anisotropic  ratios  and  constraint  factors  were  considered  for  this  data  set,  and  larger 
anisotropic  ratios  showed  some  good  improvement  of  the  data  fit  to  the  reduced  zone 
model,  indicating  that  the  low  period  data  are  more  sensitive  to  the  anisotropy  of  the 
aquifer.  But  when  considering  the  lateral  heterogeneity,  compared  to  the  pumped 
hydraulic  CPT  locations  with  surrogate  data  (26.47%  vs  11.6%  average  standard 
deviation),  the  fit  between  the  model  and  measured  data  is  not  as  good.  However,  most 
of  the  error  in  the  lateral  heterogeneity  model  occurs  around  the  six  nodes  in  the  basal 
high  K  portion  of  the  aquifer  (Fig.  79).  Fewer  multiple  intersecting  rays  through  the 
basal  zones  may  exacerbate  this;  diagonal  rays  are  needed  for  good  lateral  resolution.  If 
these  zones  are  not  considered,  the  average  percent  standard  deviation  is  lowered  to  about 
22%,  which  is  improved  and  more  comparable  to  the  error  in  the  other  pneumatic  data 
set. 

The  K  values  presented  in  Figure  80  from  HT-GP  to  HT-3  involve  some 
replacement  data  in  the  middle  and  upper  portion  of  the  test  section,  replacing  some 
unusable  field  data  (see  Lyle,  201 1  for  details).  The  ranges  of  K  values  (0.0002  to  0.001 
m/sec  [0.0006  to  0.0037ft/sec])  are  within  the  expected  range  at  GEMS.  Inversion  of  the 
experimental  phase  shift  with  an  anisotropy  ratio  of  10  produced  K  values  that  exceeded 
the  expected  range  and  several  orders  of  magnitude  greater  than  the  data  set.  A  larger 
constraint  factor  of  10  produced  reasonable  K  values  with  an  average  percent  standard 
deviation  of  about  17%;  but,  the  data  fit  tends  to  be  more  biased  towards  the  initial  K 
estimates  than  the  full  experimental  data  set  would  suggest.  In  contrast  to  the  other  CPT 
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data  sets,  the  data  fit  in  the  high  K  portion  of  the  aquifer  is  comparable  to  the  fit  in  the 
rest  of  the  aquifer. 
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Figure  79.  K  values  from  constrained  SVD  analysis  of  729  rays  at  HT-GP  (source)  to 
HT-2  (receiver).  K  values  range  from  0.0003  to  0.002  m/sec  (0.0009  to  0.0067  ft/sec). 
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Figure  80.  K  values  from  constrained  SVD  analysis  of  784  rays  at  HT-GP  (source)  to 
HT-3  (receiver).  K  values  range  from  0.0002  to  0.001 1  m/sec  (0.0006  to  0.0037  ft/sec). 
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Summary  and  Conclusions 

Since  spatial  changes  in  hydraulic  conductivity  are  a  major  factor  governing  the 
transport  and  fate  of  a  pollutant  as  it  moves  through  an  aquifer,  we  have  focused  on  the 
development  of  new  innovative  methods  to  delineate  these  spatial  changes.  The 
objective  of  the  research  presented  here  was  to  build  on  our  previous  work  to  develop  and 
improve  field  techniques  for  better  definition  of  the  three-dimensional  spatial  distribution 
of  hydraulic  conductivity  by  using  hydraulic  tomography  coupled  with  high-resolution 
slug  testing.  Hydraulic  tomography  has  received  a  lot  of  attention  in  the  literature  in 
recent  years,  but  little  field  work  has  been  done  on  the  small  scale  presented  here  using  a 
sinusoidally  varying  signal  The  first  thing  that  we  did  was  to  lay  out  the  basic  theory  of 
groundwater  flow  that  would  govern  the  experiments.  The  basic  theory  of  tidal 
influences  on  inland  wells  has  been  known  and  used  for  many  years.  The  basic  idea  is 
that  the  amplitude  and  phase  of  the  received  signal  depend  on  the  aquifer  properties 
traversed.  The  same  theory  is  applicable  to  the  research  presented  here  except  the 
geometry  and  time  scales  are  different.  Instead  of  the  ocean  providing  the  signal  over  a 
long  boundary  with  a  period  of  about  24  hours,  we  are  dealing  with  a  small  diameter  well 
providing  a  signal  over  a  relatively  short  wellbore  interval  and  having  a  period  of 
seconds.  We  developed  some  approximate  homogeneous  solutions  with  radial  geometry 
and  suggested  that  these  solutions  could  be  applied  to  the  heterogeneous  case  by  using  a 
straight  ray  spatially  weighted  approximation.  The  advantage  to  this  approach  is  that  the 
tomographic  inversion  can  be  done  without  the  need  for  an  iterative  non-linear  regression 
using  a  numerical  model.  The  traditional  approach  to  the  tomographic  inversion  is  very 
computational  intensive  and  requires  considerable  computer  resources  and  time. 
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In  order  to  process  the  data  a  number  of  new  programs  had  to  be  written.  One  of 
the  authors  (Carl  McElwee)  has  been  programming  in  Visual  Basic  for  a  number  of  years 
and  this  was  the  programming  platform  chosen.  Using  Visual  Basic  has  the  additional 
advantage  of  being  closely  interfaced  with  the  Excel  spreadsheet  program  which  is 
widely  used  and  understood.  Three  new  programs  were  written  and  improved 
continually  (denoted  by  version  numbers)  during  the  course  of  this  research.  The  first 
program  was  FitAmpPhase,  which  in  the  end  could  input  the  field  data,  take  off  a  level 
trend,  filter  the  data,  and  finally  fit  a  sine  wave  to  the  data  giving  the  amplitude  and  phase 
of  the  signal  at  the  receiver.  The  basic  theory  states  that  the  amplitude  and  phase  of  the 
received  signal  depends  on  the  diffusivity  of  the  material  through  which  the  signal 
passed.  The  diffusivity  is  the  ratio  of  hydraulic  conductivity  to  specific  storage.  After 
determining  the  amplitude  and  phase  it  was  necessary  to  extract  the  diffusivity  by  an 
inversion  process.  The  first  step  in  the  inversion  process  was  to  create  a  program  which 
would  implement  the  straight  ray  spatially  weighted  approximation.  This  program  would 
calculate  the  ray  path  lengths  in  each  zone  of  differing  K  and  calculate  the  theoretical 
phase  resulting  from  that  path.  The  output  of  this  program  (HydraulicTomAnal)  would 
be  a  set  of  matrix  equations,  with  a  matrix  of  path  lengths  multiplying  a  vector  of 
diffusivities  resulting  in  the  theoretical  phase.  An  option  for  the  HydraulicTomAnal 
program  is  to  specify  the  anisotropy  ratio  for  each  layer  The  next  step  in  the  inversion 
process  was  to  develop  a  program  that  would  use  the  matrix  of  ray  path  lengths  and  either 
theoretical  phase  or  experimental  phase  to  calculate  the  diffusivities  in  each  zone.  That 
inverse  program  (LeastSquaresSVD)  used  singular  value  decomposition  to  perform  a 
least  squares  fit  and  estimate  zonal  diffusivities.  We  usually  assumed  a  value  of  specific 
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storage  and  then  calculated  hydraulic  conductivity  from  the  diffusivity.  The  inverse 
program  could  perform  inverses  for  a  single  data  set  or  it  could  add  random  error  and 
perform  inverses  for  a  large  suite  of  data  sets  to  generate  a  Monte  Carlo  analysis. 
Another  option  for  the  inverse  program  was  to  perform  a  constrained  inverse,  where  the 
initial  input  estimates  of  K  were  given  some  additional  weight. 

Since  analytical  solutions  do  not  exist  for  the  heterogeneous  case,  numerical 
modeling  was  performed  to  check  the  spatially  weighted  straight  ray  approximation. 
Synthetic  data  were  generated  with  a  numerical  Finite  Difference  model  and  modeling 
studies  were  performed  to  demonstrate  that  data  could  be  inverted  with  reasonable 
amounts  of  error.  Zones  could  be  resolved  using  the  LeastSquaresSYD  program  to 
dimensions  of  about  1  m  (3  ft).  Comparison  of  the  4  and  30-sec  CPT  modeling  studies 
indicate  that  the  30-sec  period  CPT  does  have  more  inherent  error  than  the  3-sec  CPT, 
but  error  associated  with  the  inversion  of  field  data  is  comparable  between  the  two 
sources.  Theoretically,  it  is  expected  that  the  4-sec  data  would  be  better  than  the  30-sec 
data  due  to  averaging  over  volumes  comparable  to  the  wavelength.  The  modeling 
estimate  is  that  the  30-sec  CPT  has  about  27-29%  total  error  and  the  4-sec  CPT  has  about 
16-19%  total  error  associated  with  the  straight  ray  method,  ambient  noise,  and  the 
inversion  of  field  data. 

It  is  well  known  in  the  inverse  literature  that  the  stability  and  uniqueness  can  be 
affected  by  the  model  structure,  experimental  design,  and  data  quality.  In  our  case  this 
means  that  the  number  and  placement  of  the  zones  of  constant  K,  the  number  and 
distribution  of  ray  paths  in  our  data  sets,  and  the  ambient  noise  can  have  a  dramatic  effect 
on  the  K  results  produced  by  the  inverse  program.  We  have  performed  a  number  of 
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experimental  designs  to  arrive  at  model  structures  that  can  produce  a  reasonable  degree 
of  resolution  in  the  presence  of  expected  levels  of  noise.  In  the  end,  we  believe  that  the 
3-30  sec  period  sinusoidal  tomographic  method  applied  to  situations  like  the  GEMS  site 
is  capable  of  resolving  zones  of  K  on  the  order  of  1  m  (3  ft)  square  in  regions  having 
adequate  ray  path  coverage. 

The  application  of  the  theory  to  the  field  was  an  evolving  process;  as  the  source 
and  receiver  equipment  was  refined,  the  field  procedures  needed  to  be  reassessed.  The 
HRST  procedures  and  equipment  used  by  other  researchers  at  the  GEMS  location  guided 
the  transition  to  sinusoidal  CPT  and  made  it  relatively  easier.  Many  of  the  equipment 
items  used  in  CPT  were  initially  borrowed  or  adapted  from  previously  existing  HRST 
equipment.  To  make  this  type  of  testing  readily  available  to  future  researchers,  the 
majority  of  the  equipment  was  put  together  from  “off  the  shelf’  sources.  After  an 
experiment  was  conducted,  the  equipment  could  be  revamped  to  suit  the  particular 
stratagem  of  future  experiments.  It  was  important  to  us  to  design  equipment  that  could  be 
used  in  standard  2  inch  monitoring  wells,  so  compact  design  was  important.  The  single 
source  location  was  isolated  by  a  straddle  packer  inflated  with  nitrogen.  Receiver 
location  varied  from  one  to  five  for  a  given  record,  depending  on  whether  we  were  doing 
a  ZOP  or  MOG  well  survey.  In  any  case,  each  receiver  location  was  also  isolated  by 
straddle  packers.  Early  research  results  showed  the  importance  of  packing  off  the 
measurement  interval  to  obtain  the  best  quality  data.  In  the  course  of  this  research 
several  designs  for  multilevel  receivers  and  types  of  pressure  transducers  were  evaluated, 
with  several  resulting  failures.  In  the  end,  we  developed  two  multilevel  receivers  each 
containing  five  pressure  transducers  that  yielded  good  data.  The  source  and  receiver 
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locations  were  established  in  absolute  elevation  by  the  use  of  steel  tapes  attached  to  the 
riser  pipes  measuring  the  distance  below  top  of  casing,  and  then  surveying  in  all  the  well 
tops  with  respect  to  a  benchmark  at  the  site.  Two  mechanisms  were  used  to  obtain 
oscillatory  source  inputs.  The  first  used  air  pressure  from  a  compressor  to  depress  and 
oscillate  the  water  column  in  the  source  well.  It  was  necessary  to  stay  close  to  the  natural 
frequency  of  the  well  to  obtain  a  good  sinusoidal  signal,  so  a  period  in  the  3-4  sec  range 
was  used.  The  second  source  used  mechanical  pumping  with  a  computer  controlled 
servo-valve  to  inject  water  into  the  aquifer  to  create  a  sinusoidal  pressure  signal.  We 
were  limited  by  available  equipment  to  a  30  sec  period  for  this  source.  Additionally,  an 
in-situ  pneumatic  source  with  a  computer  controlled  3-sec  period  was  developed  and 
deployed  by  a  Geoprobe  rig  to  evaluate  the  feasibility  of  tomographic  application  by 
direct-push  technology.  The  two  multi-level  receivers  were  used  in  separate  wells  so  data 
acquisition  was  very  efficient  while  using  the  Geoprobe  source. 

Six  new  wells  were  installed  at  GEMS  for  this  project.  The  wells  were  installed 
with  direct  push  technology  (Geoprobe),  which  causes  less  aquifer  disturbance  than  many 
other  methods.  Electrical  conductance  profiles  (ECPs)  were  obtained  at  the  first  three 
well  locations  and  we  were  able  to  discern  the  boundaries  between  the  Tonganoxie 
Sandstone  member,  the  unconsolidated  sands  and  gravels,  and  the  confining  silts  and 
clays.  The  locations  were  chosen  to  give  an  area  of  study  with  some  lateral  extent. 
Basically  there  is  a  center  well  and  five  surrounding  wells  allowing  for  five  cross  sections 
connecting  the  center  well  and  surrounding  ones.  The  22  m  (70  ft)  deep  wells  were  fitted 
with  .05  m  (2  in)  PVC  casing  and  with  11m  (35  ft)  of  screen  at  the  bottom.  These  wells 
were  extensively  developed  to  prepare  them  for  data  collection.  Three  of  the  existing 
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wells  at  GEMS  were  also  used  for  early  ZOP  work.  High  resolution  slug  testing  (HRST) 
was  performed  on  each  of  the  wells  used  in  this  research  at  .30  m  (1  foot)  intervals.  The 
HRST  vertical  profiles  of  K  values  at  each  well  were  invaluable  for  delineating  high  and 
low  K  zones  and  constraining  the  tomographic  inverses. 

Early  data  collection  on  the  wells  consisted  of  ZOP  surveys  looking  only  at 
horizontal  ray  paths;  but  later  work  consisted  of  MOG  surveys  with  many  diagonal  ray 
paths.  In  the  beginning  two  types  of  source  geometries  were  used:  the  whole  well  line 
source  (entire  open  screen  water  column  oscillated)  and  the  isolated  point  source 
(oscillation  was  only  in  contact  with  a  small  packed  off  interval  of  the  screen).  The  line 
source  introduces  a  greater  amount  of  energy,  and  therefore  has  a  greater  propagation 
distance  than  the  point  source.  On  the  other  hand,  the  point  source  geometry  allows  for  a 
better  vertical  resolution  of  the  aquifer  characteristics.  For  the  ZOP  work  seven 
pneumatic  line  source  CPT  profiles  and  six  point  source  CPT  profiles  were  completed  at 
GEMS.  MOG  surveys  were  completed  between  all  five  well  pairs  each  involving  the 
center  well,  both  using  pneumatic  3  to  4-sec  data  and  mechanical  pumping  30-sec  data. 

In  addition,  we  adapted  the  pneumatic  3 -sec  source  method  to  be  used  with  a  Geoprobe 
unit.  We  collected  MOG  data  in  two  receiver  wells  simultaneously  while  the  Geoprobe 
unit  pulled  the  source  up  in  one  foot  increments  from  the  bottom  of  the  aquifer.  The 
range  of  radial  distances  between  tested  well  pairs  was  1.5  to  1 1.5  m. 

The  horizontal  ray  (ZOP)  data  that  were  collected  show  that  the  CPT  K  profiles 
mimic  the  general  trends  in  the  HRST  K  profiles  measured  at  the  respective  wells. 
Overall,  the  CPT  data  appear  to  average  the  K  profiles  of  the  well  pair  in  question. 
However,  there  are  important  differences.  The  heterogeneities  of  the  aquifer  between  the 
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well  pair  are  probably  the  cause  of  this  difference,  and  the  difference  can  not  be  fully 
explained  without  collecting  diagonal  ray  path  data,  using  more  advanced  models,  and 
using  tomographic  methods.  We  made  an  elementary  attempt  at  calculating  “anomalous 
K  values”  between  well  pairs.  According  to  simple  theory,  the  CPT  data  should  be  an 
average  of  the  HRST  measurements  and  should  fall  somewhere  within  the  HRST  limits. 

If  not,  then  the  differences  may  represent  heterogeneities  between  the  well  pair.  Both  the 
line  source  data  and  the  point  source  data  appear  to  distinguish  variations  in  K  that  are 
not  present  in  the  HRST  data.  However,  the  point  source  data  appear  to  have  the  best 
resolution  of  the  data  presented  here.  The  line  source  method  is  very  difficult  to  interpret 
and  was  not  used  in  later  experiments.  Experiments  indicate  that  the  data  are 
reproducible  at  different  times  and  with  the  source  and  receiver  wells  reversed,  within 
experimental  accuracy.  Initially,  both  amplitude  and  phase  data  were  used  to  estimate 
CPT  K  values;  however,  amplitude  data  are  less  reliable  than  the  phase  data  because  of 
the  exacting  experimental  measurements  that  are  necessary.  Therefore,  relative 
amplitude  data  were  only  used  as  a  diagnostic  tool  and  not  used  for  determining  K  in  later 
work.  The  ZOP  CPT  method  can  not  estimate  lateral  heterogeneity,  but  it  can  be  of 
practical  use  for  discerning  changing  horizontal  flow  units. 

Five  well  pairs  were  analyzed  using  pneumatic  CPT  MOG  data  with  a  period  of  3 
to  4-sec.  The  MOG  data  sets  allowed  the  full  use  of  tomographic  methods  to  determine 
the  K  distribution.  All  had  reasonable  interwell  K  distributions  after  using  a  constrained 
inversion,  as  compared  to  the  general  range  seen  with  HRST.  K  values  at  the  site  are 
known  from  HRST  to  range  from  about  0.000305  m/s  to  0.00305  m/s,  and  follow  the 
general  trend  of  a  higher  K  zone  at  the  base  of  the  aquifer,  a  low  K  zone  above,  a 
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moderately  high  K  zone  still  higher  up  in  the  profile,  and  a  low  K  zone  at  the  top.  The 
success  or  failure  of  the  inversion  was  evaluated  by  comparing  the  resulting  K  values  to 
the  range  of  K  values  seen  from  HRST  as  well  as  the  general  trends  of  high  or  low  K 
zones  seen  from  HRST.  The  success  of  the  inversion  seems  to  be  correlated  with  the 
number  of  ray  paths  between  the  source  and  receiver  wells.  Varying  source  and  receiver 
intervals  for  each  well  pair  offered  the  opportunity  to  examine  how  much  data  needed  to 
be  collected.  Initial  data  were  collected  at  a  fine  scale  (0.305  m)  given  the  resolution 
capabilities  of  the  model.  Editing  the  number  of  ray  paths  used  in  the  inversion  process 
allowed  comparisons.  Of  the  variations  tested  in  this  study,  the  geometry  used  for  GEMS 
was  most  efficiently  and  accurately  characterized  with  about  300  ray  paths,  but  750  ray 
paths  will  provide  some  additional  accuracy  if  time  is  available  for  their  collection. 

Some  small  problems  with  some  well  pairs  can  be  explained  by  equipment  problems  in 
one  case  and  by  too  few  ray  paths  in  another  two  cases.  The  results  of  this  research  show 
that  hydraulic  tomography  combined  with  appropriate  inversion  programs  can  estimate 
interwell  K  distributions  with  resolutions  down  to  about  one  square  meter  in  the  most 
sensitive  regions. 

Pumped  hydraulic  CPT  MOG  data  with  a  30-sec  period  were  collected  from  a 
radial  well  array  with  a  central  source  well  (HT-3)  and  five  receiver  wells  (HT-1,  HT-2, 
HT-4,  HT-5  and  HT-6).  In  addition,  the  pneumatic  method  with  a  period  of  3-sec  was 
used  with  the  Geoprobe  as  source  to  collect  MOG  data  in  wells  HT-2  and  HT-3.  This 
generated  a  large  tomography  data  set  which  was  used  to  evaluate  the  variation  of  K 
across  the  area.  In  some  instances,  poor  data  were  collected  due  to  equipment  failure  or 
poor  signal  to  noise  ratio.  To  address  that  data  loss,  a  method  was  devised  to  replace  the 
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bad  data  with  data  based  on  the  HRST  K  point  source  data,  so  that  the  larger  data  sets  of 
multiple  MOGs  in  a  well  pair  could  be  preserved  and  tomographically  analyzed  as  a 
whole.  Various  models  of  isotropy,  anisotropy,  and  lateral  heterogeneity  were  used  to 
evaluate  the  goodness  of  the  fit  coming  from  the  inversion  process.  Different  degrees  of 
constraint  were  evaluated  to  determine  the  best  data  fit  by  SVD  analysis.  Data  analysis 
indicated  that  the  3 -sec  period  data  were  more  sensitive  to  vertical  anisotropy  than  had 
been  expected.  Modeling  with  a  greater  degree  of  anisotropy  allowed  a  better  fit  for  K 
values  using  the  3 -sec  data  when  considering  fine-grained  material  and  long-offset  rays, 
while  the  30-sec  CPT  period  data  seemed  somewhat  insensitive  to  different  degrees  of 
anisotropy.  In  general  the  3-sec  data  needed  more  constraint  to  produce  a  stable  inverse. 
The  ranges  of  calculated  K  values  were  compared  to  the  lithology  of  the  aquifer  and 
HRST  K  values  from  the  well  array  to  evaluate  the  success  of  the  inversion. 

Tomographic  analysis  from  this  research  generated  K  values  that  fell  within  these 
guidelines,  indicating  good  performance  of  the  CPT  equipment  and  data  processing 
techniques.  The  hydraulic  pumping  30-sec  CPT  data  were  compared  to  previous  research 
completed  with  a  pneumatic  3-4  sec  CPT  period.  Data  trends  and  K  values  were  similar 
and  within  the  general  range  seen  with  HRST,  although  the  30-sec  data  did  not  have 
some  of  the  resolution  obtained  with  the  shorter  period  CPT  source  used  for  earlier 
research  data  and  the  direct  push  data. 

This  project  has  provided  support  for  three  students  to  finish  the  Masters  Degree 
and  some  support  for  four  others  during  their  undergraduate  or  graduate  degrees.  This 
work  has  shown  that  the  use  of  an  oscillatory  source  together  with  tomographic 
techniques  and  supporting  high  resolution  slug  testing  is  a  promising  hydrogeological 
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tool  to  determine  interwell  hydraulic  conductivity  distributions.  The  data  sets  obtained  in 
the  course  of  this  project  have  provided  hydraulic  conductivity  estimation  in  a  360-degree 
radial  array  over  an  extended  area  at  GEMS.  We  have  shown  that,  at  least  at  GEMS,  a 
resolution  of  about  one  meter  square  can  be  achieved  in  areas  with  adequate  ray  path 
coverage.  We  have  designed  equipment  that  allows  this  work  to  be  done  in  standard  2 
inch  monitoring  wells  with  long  screens.  We  have  also  shown  that  the  work  can  be  done 
using  a  direct  push  unit  (Geoprobe)  to  deploy  the  source.  We  were  constrained  to  using 
signals  with  periods  of  about  3-sec  and  30-sec  due  to  equipment  limitations.  We  know 
that  the  resolution  is  dependent  upon  period,  so  further  research  needs  to  be  done  using  a 
pumping  source  that  can  bridge  the  gap  between  3  and  30  seconds.  This  research  was 
supported  in  part  by  the  U.S.  Department  of  Defense,  through  the  Strategic 
Environmental  Research  and  Development  Program  (SERDP). 
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Appendix  B:  HRST  K  Profiles 


Well  HT-1  HRST 


K(itYs) 


Figure  Bl:  HRST  results  for  well  HT-1  (processed  by  Brett  Engard). 
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Well  HT-2  HRST 


Figure  B2:  HRST  results  for  well  HT-2  (processed  by  Brett  Engard). 
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Figure  B3:  HRST  results  for  well  HT-3  (processed  by  Brett  Engard). 
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VNfell  HT-4  HRST 


Figure  B4:  HRST  results  for  well  HT-4  (processed  by  Pema  Deki). 
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Well  HT-5  HRST 


Figure  B5:  HRST  results  for  well  HT-5  (processed  by  Pema  Deki). 
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Well  HT-6  HRST 


Figure  B6:  HRST  results  for  well  HT-6  (processed  by  Pema  Deki). 
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